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^HOOL 
I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  problem  of  motion  stability  of  ships  and  other  marine  vehicles  has  been 
the  subject  of  extensive  studies  in  the  past  (Comstock,  1977).  Most  of  these 
studies  are  with  regards  the  ship's  directional  stability  in  open  waters  under 
open  loop  conditions.  It  is  well  known  that  in  restricted  waters  such  as  canals 
or  rivers,  although  it  is  possible  to  have  positional  stability  in  principle,  in 
reality  it  is  never  the  case.  This  is  due  to  the  destabilizing  effects  of  the  bank 
suction  forces  and  moments.  These  act  always  in  a  destabilizing  fashion;  i.e., 
after  a  small  initial  disturbance  they  produce  a  force  or  a  moment  which  tends 
to  increase  the  ship's  path  deviation  from  nominal  (Comstock,  1977).  Open 
loop  conditions,  important  as  they  are,  rarely  represent  real  life  applications. 
Ships  traveling  in  canals  are  under  closed  loop  control,  typically  provided 
through  the  helmsman  or  an  autopilot.  It  becomes  then  important  to  assess 
the  stability  of  the  system  under  finite  disturbances  and  incorporating  some 
modeling  of  closed  loop  operations. 


B.  OBJECTIVES  AND  OUTLINE 

This  thesis  utilizes  both  linear  and  nonlinear  techniques  in  order  to  analyze 
the  problem  of  closed  loop  dynamic  stability  of  ships  in  canals.  Ship  modeling 


is  discussed  in  Chapter  II  and  it  follows  standard  ship  maneuvering  equations 
(Clayton  and  Bishop,  1982)  augmented  by  a  model  for  bank  suction  effects 
(Beck,  1976).  A  Mariner  class  ship  is  selected  as  a  model  because  of  the  avail- 
ability of  data  for  its  hydrodynamic  coefficients  (Comstock,  1977),  (Bernitsas 
and  Kekridis,  1984).  A  heading  control  law  based  on  Nomoto's  first-order 
model,  coupled  with  a  line  of  sight  guidance  law  is  chosen  in  order  to  model 
the  fundamental  behavior  of  closed  loop  conditions.  Since,  in  practice  control 
actions  are  hardly  ever  instantaneous,  appropriate  time  lags  are  introduced 
in  the  feedback.  These  are  modeled  using  the  techniques  outlined  in  (Venne, 
1992).  Linear  eigenvalue  analysis  is  performed  in  Chapter  III  in  order  to 
reveal  parametric  stability  boundaries.  It  is  established  that  loss  of  stability 
occurs  when  a  pair  of  complex  conjugate  eigenvalues  crosses  the  imaginary  axis 
which  results  in  periodic  solutions  or  limit  cycles  (Guckenheimer  and  Holmes, 
1983).  A  nonlinear  analysis  based  on  Hopf  bifurcation  methods  (Chow  and 
Mallet-Paret,  1977)  and  (Hassard  and  Wan,  1978)  is  performed  in  Chapter 
IV.  Conclusions  derived  from  this  study  and  some  recommendations  for  fur- 
ther extensions  are  discussed  in  Chapter  V.  The  basic  programs  that  were  used 
to  produce  the  results  are  included  in  the  Appendix. 


II.  PROBLEM  FORMULATION 


A.  SHIP  DYNAMICS 

Restricting  our  attention  to  the  horizontal  plane,  the  mathematical  model 
consists  of  the  nonlinear  sway  (translational  motion  parallel  to  the  vehicle's 
longitudinal  axis)  and  yaw  (rotational  motion  about  the  vertical  axis  )  equa- 
tions of  motion.  If  we  consider  a  moving  Cartesian  coordinate  frame  located 
at  the  geometric  center  of  the  vehicle,  Newton's  equations  of  motion  are: 

m(v  +  ur  +  xcf)    =    Y,  (1) 

Izt  -\-mxaiv  +  ur)    =    N  ,  (2) 

where  v  and  r  are  the  relative  sway  and  yaw  velocities  of  the  moving  vehicle 
with  respect  to  the  water,  m  is  the  vehicle's  mass,  Iz  is  its  moment  of  in- 
ertia with  respect  to  the  vertical  axis,  and  u  is  the  constant  forward  speed. 
Since  all  quantities  will  be  considered  as  dimensionless  in  this  study,  in  the 
following  we  consider  u  to  be  equal  to  one.  xq  is  the  longitudinal  position 
of  the  ship's  center  of  gravity  with  respect  to  its  centroid,  and  Y,N  represent 
the  total  excitation  sway  force  and  yaw  moment  respectively.  Following  stan- 
dard ship  maneuvering  assumptions,  these  forces  can  be  expressed  as  the  sum 
of  quadratic  drag  terms  and  first-order  memoryless  polynomials  in  v  and  r, 
which  represent  added  mass  and  damping  due  to  water.  In  this  way  they  can 


be  expressed  by: 

Y    =    YyV  +  Yrf  +  Yvv  +  -Yvvvv3  +  -Yvrrvr2  + 

-Yrvvrv2  +  Yrr  +  Yfy,  y)  +  YSS  ,  (3) 

1  1 

N    =    NyV  +  NfT  +  Nvv  +  -Nvvvv3  +  -Nvrrvr2  + 

-Nrvvrv2  +  Nrr  +  N{^,y)  +  N66  ,  (4) 

where  Ya,  Na  represent  the  partial  derivatives  with  respect  to  the  indicated 
variable  a  and  8  is  the  effective  rudder  angle.  Y(ip,y),  N(ip,y)  represent  the 
interaction/proximity  forces  and  moments  that  arise  due  to  the  presence  of 
the  canal,  and  shallow  water  effects. 

The  resulting  nonlinear  differential  equations  can  be  non-dimensionalized 
with  respect  to  the  constant  forward  speed  of  the  ship,  u  and  its  length,  I. 
The  dimensionless  time  variable  is  then  equal  to  tu/l.  A  standard  inertial  sys- 
tem (x,y)  is  introduced  here  where  the  a: -axis  points  in  the  assumed  nominal 
straight  line  path  and  the  t/-axis  is  the  distance  from  this  nominal  path,  see 
Figure  1.  We  assume  that  the  nominal  straight  line  path  is  along  the  center- 
line  of  the  canal.  In  cases  where  the  concept  of  a  geometric  centerline  is  not 
applicable,  we  assume  that  the  nominal  path  is  along  the  zero  bank  suction 
location.  This  is  consistent  with  recommended  navigation  practices  that  are 
currently  in  use. 

The  complete  model  is  presented  by  the  following  equations: 

(m  -  Yy)v  -  (Yr  -  mxG)r  =  Yvv  +  -Yvvvv3  +  -Yvrrvr2  + 


2^ 


Canal  centerline       <- 


Figure  1:  Vehicle  geometry  and  definitions  of  symbols 


-Yrvvrv2  +  (Yr  -  m)r  +  Yty,  y)  +  Y66  , 

-(Ny  -  mxG)v  +  (Iz  -  Nr)r  =  Nvv  +  -Nvvvv3  +  -Nvrrvr2+  , 

b  2 


1 


Nrvvrv2  +  {Nr  -  mxG)r  +  N(ip,  y)  +  NSS  , 


2 
ip  =  r  , 

y  =  sin  tp  +  v  cos  ijj  , 


(5) 

(6) 

(?) 
(8) 


where  the  expressions  for  the  ship's  yaw  rate  as  well  as  its  inertial  position  rate 
have  been  added,  and  tp  is  the  local  heading  angle.  The  interaction/proximity 
forces  and  moments  are  expanded  to  include  up  to  third  order  terms, 


1  3       1  3 

Y (^,  y)  =  Y^tp  +  Yyy  +  -Y^ip    +  -Yyyyy    , 

N{if>,  y)  =  N^  +  Nyy  +  -N^^3  +  -Nyyyy3  , 

o  o 


(9) 
(10) 


as  explained  in  more  detail  in  the  folowing  section. 


B.  HYDRODYNAMIC  COEFFICIENTS 

We  chose  a  Mariner  class  ship  as  a  representative  model.  Its  hydrodynamic 
coefficients  and  geometric  properties  were  taken  from  (Comstock,  1977)  and 
(Bernitsas  and  Kekridis,  1984).  Results  from  (Beck,  1976)  were  used  in  order 
to  model  the  bank  suction  forces  and  moments.  Typical  results  are  shown  in 
Figure  2.  These  show  force  and  moment  coefficients  versus  ship  deviation  (77), 
and  for  canal  width  w  =  OAL.  Force  and  moment  coefficients  are  nondimen- 
sionalized  with  respect  to  the  water  density  and  the  ship's  speed  and  length, 
as  is  customary  in  ship  maneuvering.  Since  the  suction  force  and  moment 
must  change  their  sign  as  either  y  or  ij;  changes  its  sign,  they  must  be  modeled 
by  odd  power  polynomials,  as  was  done  in  the  previous  section.  Numerical 
values  for  the  coefficients  Yy,  Yyyy,  Ny,  and  Nyyy  can  be  found  by  curve  fitting 
of  the  results  of  Figure  2.  Using  a  depth  to  draft  ratio  of  1.9  we  were  able  to 
estimate, 

Yy    =    0.02  , 
Yyyy    =    0.468, 
Ny    =    -0.0025  , 

Nyyy      —      0  . 

The  value  of  N^  was  estimated  by  "discretizing"  the  ship  in  two  segments, 


FIGURE    13 

VARIATION   OF    SIDE    FORCE    AND    MOMENT    WITH 

WALL    POSITION    RATIO    FOR   THE   MARINER 


"      -1, 


-2.  - 


-3 

-4. 


—  c      _ 


-6. 


w/L   =    .4 


U   =   7   kts,    full    scale 

F.     =    .37,  h/T   =   1.3 

n  ' 

=    .31, h/T   =1.9 


Experiment     


Figure  2:  Forces  and  moments  due  to  canal  [Beck  (1976)] 

7 


at  the  bow  and  the  stern,  and  calculating  the  suction  forces  on  each  part 
individually.  Their  resultant  moment  produced, 

Nj,  =  0.01  . 

The  value  of  N^^  was  taken  to  be  zero,  due  to  lack  of  reliable  data.  The  value 
of  Y$  was  estimated  to  be  equal  to  —  Yv  =  0.014  which  is  true  for  motions  along 
the  centerline  of  the  canal  (Comstock,  1977).  Again,  due  to  lack  of  reliable 
data  we  took  Y^^  =  0. 


C.  CONTROL  LAW 

The  linearized  set  of  equations  (5)  through  (7)  can  be  expressed  in  the 
following  form: 

ip    =    r  , 

v     =    anip  +  a\iv  +  a^r  +  auy  +  b\6  ,  (11) 

r     =    a2i?/>  +  a22v  +  a23r  +  a24y  +  b26  , 

where  the  coefficients  a^,  bi  are  functions  of  the  vehicle  hydrodynamic  coeffi- 
cients and  geometric  properties  and  are  given  below: 

an    =     —  [(I z  -  Nr)Yj,  +  (Yf  -  mxG)N^}  , 

Uyr 

al2    =     —  [{Iz  ~  Nr)Yv  +  (Y+  -  mxG)Nv]  , 

ai3    =     — —  [{Iz-  Nr)(Yr  -  m)  +  (Y+  -  mxG){Nr  -  mxG))  , 

L> vr 

a14     =     —  [{Iz-Nr)Yy  +  (Yr-mxG)Ny\  , 


Q>21 
0-22 

«24 
&1 

b2 

D  1)T 


D, 


D, 


Dv 

1 


u vr 
1 

L) Vr 
1 


[(Nv  -  mxG)Y^  +  (m  -  iyj\fy]  , 

[(TV;,  -  mic)^  +  (m  -  3^)iVj  , 

[(No  -  mxcOO'r  -  m)  +  (m  -  Yv)(Nr  -  mxG)]  , 

[(Nv  -  mxG)Yy  +  {m-  Yv)Ny]  , 

[{Iz  -  N,)Y6  +  (Y,  -  mxG)N6]  , 

[(Nv  -  mxG)Y6  +  (m  -  YV)N6]  , 


(m  -  YV){IZ  -  Nr)  -  {Nv  -  mxG)(Yr  -  mxG) 


A  control  law  that  could  model  a  human  operator  should  not  be  based  on 
feedback  of  the  side  slip  velocity  v.  Instead,  it  is  more  likely  that  human 
operators  will  respond  to  changes  in  the  ship's  heading  angle  ip  and  rate  of 
change  of  heading,  r.  Therefore,  we  choose  to  base  our  control  law  on  Nomoto's 
model,  which  follows  by  assuming  negligible  sway  velocity  v , 


r  =  ar  +  cip  +  b6 


(12) 


where  the  coefficients  are 


a    =    a23  , 
b    =    b2, 
c    =    a2\  ■ 

A  linear  control  law  can  be  introduced  in  the  form, 


60  =  ki{ip  -  ipc)  +  k2r  , 


(13) 


where  ipc  is  the  commanded  heading  angle  and  the  gains  fci,  k2  are  computed 
such  that  the  closed  loop  system  (7),  (12),  and  (13)  has  the  desired  dynamics. 
The  existence  of  the  difference  (ip  —  ipc)  in  the  control  law  (13)  has  the  effect  of 
trying  to  point  the  ship's  longitudinal  axis  towards  the  desired  heading.  The 
characteristic  equation  of  the  system  is  obtained  from  (7),  (12),  and  (13)  as 

s2  -  {a  +  bk2)s  -  (c  +  bki)  =  0  . 

If  the  desired  characteristic  equation  is 

S2  +  2(LUnS  +  uj2n  =  0  , 

the  control  gains  can  be  computed  from 

6        6' 

a  +  2Co>n 

ko    = • 

b 

The  natural  frequency  ton  and  damping  ratio  £  are  selected  based  on  general 

properties  of  second  order  systems.  Relatively  high  values  of  u>n  and  low  values 

of  £  will  correspond  to  a  responsive  operator,  while  the  opposite  is  true  for 

combinations  of  low  un  and  high  (. 

Finally,  in  order  to  take  into  account  the  effect  of  rudder  saturation,  the 
commanded  rudder  angle  is  given  by 

<5  =  <5sattanh(^   ,  (14) 

where  <5o  is  the  slope  of  6  at  the  origin  given  by  (13),  and  <5sat  is  the  saturation 
limit  on  6,  typically  around  0.4  radians.  The  hyperbolic  tangent  function  is 
used  instead  of  a  hard  saturation  function,  because  of  its  differentiability. 
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D.  GUIDANCE  SCHEME 

Since  the  previous  control  law  stabilizes  the  ship  to  any  commanded  head- 
ing angle,  it  must  be  coupled  with  an  appropriate  orientation  guidance  scheme 
to  provide  path  keeping  along  the  desired  track.  The  simplest  guidance  scheme 
which  models  some  fundamental  aspects  of  helmsman  behavior  is  a  pure  pur- 
suit guidance  where  the  commanded  heading  angle  equals  the  line  of  sight 
angle, 

i/>e=-  tan"1  (—)  ,  (15) 

\Xd/ 

as  shown  in  Figure  1.  The  ship  is  located  at  (x,y)  and  attempts  to  point  its 
longitudinal  axis  towards  a  target  point  which  is  located  ahead  of  the  ship 
on  the  reference  path  at  a  constant  preview  distance  Xd-  Pursuit  guidance  is 
achieved  by  commanding  a  heading  angle  ipc  equal  to  the  line  of  sight  angle 
(15). 

The  positional  error  information  y  in  equation  (15),  is  assumed  to  lag  the 
actual  error  y  by  an  amount  of  T  seconds,  in  other  words, 

y  =  y(t-T).  (16) 

In  this  equation,  the  time  lag  T  models  the  necessary  time  that  it  takes  for 
the  helmsman  to  process  his  path  deviation  and  initiate  appropriate  corrective 
actions. 
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III.  LINEAR  ANALYSIS 

In  this  Chapter,  we  present  a  linearized  analysis  of  the  equations  of  motion. 
The  purpose  of  this  analysis  is  to  assess  stability  or  instability  of  the  equations 
in  response  to  small  deviations  from  the  straight  line  reference  track.  No 
attempt  will  be  made  here  to  characterize  the  transient  response  of  the  system. 
This  can  be  accomplished  through  numerical  simulations. 


A.  COMBINED  SYSTEM 

If  we  incorporate  the  interaction/proximity  forces  (9)  and  (10)  in  the  ship 

dynamics  model,  equations  (5)  through  (8),  we  end  up  with  the  combined 

system, 

tp  =  r  , 

(m  -  Yy)v  -  {Yf  -  mxG)r  =  Yvv  +  -Yvvvv3  +  -Yvrrvr2  +  -Yrvvrv 2+ 

1  3       1  3 

(Yr  -  m)r  +  Y^tp  +  Yyy  +  -Y^V    +  ^YyyyV    +  W  , 

1  1 

-{Ny  -  mxG)v  +  (Iz  -  Nr)r  =  Nvv  +  -Nvvvv3  +  -Nvrrvr2+  (17) 

1  6  1  2 

-Nrvvrv2  +  (Nr  -  mxc)r  +  N^ip  +  Nyy  +  -N^^ip3+ 


-Nyyyy3  +  Ns6 


I3 

y  =  sin  tp  +  v  cos  i/j  . 

Study  of  the  asymptotic  properties  of  this  system  is  the  subject  of  this  and 
the  following  chapters. 
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B.  LINEARIZATION 

The  linearized  form  of  equations  (17)  is  the  following, 

ip  =  r  , 

(m  —  Yy)v  —  (Yj-  —  mxo)r  =  Yvv  +  (Yr  —  m)r+ 

Y^  +  Yyy  +  Ys6,  (l8) 

-(Nv  -  mxG)v  +  (Iz  -  Nr)r  =  Nvv  +  (Nr  -  mxG)r+  l     ' 

N^ip  +  Nyy  +  N66  , 
y  =  ip  +  v  . 

The  rudder  angle  8  has  one  of  the  forms  that  are  developed  below.  The  time 
delay  operator  can  be  expressed  in  terms  of  its  Taylor  series  expansion, 

y(t  -  T)  =  y  -  Ty  +  l-T2y  -  ^T3y^  +  ■■■  .  (19) 

Practical  computations  can  be  performed  by  truncating  equation  (19)  to  first, 
second,  or  third  order. 

1.  First  order  approximation  in  y 
We  have, 

y(t  -  T)  =  y  -  Ty  , 


or 


In  its  linear  form, 


and,  therefore, 


y(t-T)=y-T(iP  +  v).  (20) 

tanh  (  T2" )  =  T2"  ' 


6  =  80. 
Furthermore,  in  its  linear  form  equation  (15)  can  be  written  as 

Ipc  = • 
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If  we  incorporate  (20)  into  (13),  we  derive  the  linearized  first  order  approxi- 
mation in  8, 

8  =  kiip  +  k2r  +  —  y (ip  +  v)  .  (21) 

Xd  xd 

2.  Second  order  approximation  in  y 

Keeping  the  second  order  terms  in  y(t  —  T)  we  get, 

y{t-T)=y-Ty  +  ^T2y  .  (22) 

If  we  incorporate  (22)  into  (13)  along  with  the  linear  equations  8  =  <5o  and 
ipc  =  -y(t  -  T)/xd,  we  get 

8  =  kill)  +  k2r  +  —  y (ip  +  v)  +  -^ — (r  +  v).  (23) 

Xd  Xd  2xd 

3.  Third  order  approximation  in  y 


In  this  case, 


y(t-T)=y-Ty  +  ±T2y  -  ^Tzy^  ,  (24) 


and 


h         hT  .  k{T2  kiT3  .  .  .     . 

8  =  kiip  +  k2r  H y (ip  +  v)  +  — — (r  +  v) (r  4-  v)  .        (25) 

Xd  Xd  2xd  oxd 


4.  First  order  approximation  in  8 

If  a  time  lag,  T,  exists  on  8  instead  of  simply  y,  then 

8  =  <5sat  tanh  f  — !-  J   , 

where, 

8i  =  80{t  -T)=80-  T80  . 
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This  equation  models  a  time  lag  associated  with  the  entire  application  of  the 
necessary  corrective  control  action  and  not  just  its  positional  error.  Using  the 
above  equations  along  with  the  linear  equations  8  =  8\  and  ipc  =  —y/xd,  we 

get 

k\  k{T 

8  —  kiip  +  k2r  H y  —  k\Tr [ip  +  v )  —  k2Tr  .  (26) 

Xd  Xd 

5.  First  order  approximation  in  both  y  and  8 

In  a  more  general  setting,  we  can  assume  that  a  time  lag,  7\,  exists  in  the 
control  law  8,  and  a  different  time  lag,  Tb,  is  present  in  the  processing  of  the 
positional  error  y.  Assuming  a  first  order  approximation  for  both,  we  have 

<5i  =  8o(t  -  Ti)  =  6o  -  Ti80  , 

and 

y(t  -T2)  =y  -  T2y  . 

Therefore,  in  this  case  using  the  above  equations  along  with  the  linear  equa- 
tions 8  =  8\  and  ipc  =  —y{t  —  l2)/£d  ,we  obtain 

k\  k{T\ 

8  =  kiip  +  k2r  -\ y  —  k{Tir (ip  +  v)  — 

Xd  xd 

bTl^  +  ^  +  bMlir  +  i,)-^.  (27) 

Xd  Xd 

Using  one  of  the  above  expressions,  the  linearized  equations  of  motion  can 

be  written  as 

ip  =  r  , 

v  =  anip  +  ai2v  +  aur  +  o14y  +  biS  ,  ^g) 

r  =  a2iip  +  a22v  +  a2zr  +  a24y  +  b28  , 

y  =  ip  +v  , 

where  all  coefficients  have  been  previously  defined. 
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C.  SERIES  EXPANSIONS  OF  TIME  LAG 


1.  First  order  approximation  in  y 


In  this  case  we  have, 


y(t-T)=y-Ty. 


Substitution  of  (21)  into  (28)  yields  the  following  matrix  system, 


V 

r 

y  J 


0  0  10 

-421,1      -422,1      -423,1      -424,1 
-431,1      -432,1      ^33,1      -434,1 

110  0 


V 

r 

y 


(29) 


where, 


-421,1 

= 

an  +  b\ki  — 

-422,1 

= 

an 

Xd 

-423,1 

= 

ai3  +  bik2 

-4  24,1 

= 

biki 

Oi4  + 

Xd 

-431,1 

= 

021  +  &2&1  _ 

-432,1 

= 

b2k{T 

0-22 

Xd 

-433,1 

= 

£23  +  ^2^2 

-434,1 

= 

b2h 

a24  + 

bikxT 

Xd 


Xd 


Xd 


2.  Second  order  approximation  in  y 


In  this  case  we  have, 


2-: 


y(t-T)=y-Ty  +  -Tly 
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Substitution  of  (23)  into  (28)  yields  the  following  matrix  system, 


[ip] 

V 



r 

.  y  . 

where, 


-421,2  = 
^22,2  = 
^23,2      = 


0  0  10 

^21,2      ^22,2      ^23,2      ^24,2 
^31,2     ^32,2      ^33,2      ^34,2 

110  0 


auXd  +  hkiXd  —  bikiT 

xd  -  0.5&i/ciT2 
ai2Xd  ~  b\k\T 
xd  -  O.Sb^T2 
auXd  +  bik2xd  +  0.5bjkiT2 


4> 
v 
r 

y 


xd  -  O.bbfaT2 
auxd  +  b\k\ 

24,2    ~     xd  -  0.56ifciT2 

b2k1T       b2k1T2 

Mia   —   a2i  +  &2«i  - 


xd  2xd 

>2 


^21,2 


^32,2      —      0-22 


b2k{T       b2kiT 

■  —Z ^22,2 


Xd 


•^33,2      =      0-2Z  +  b2k2  + 


2xd 
b2kiT2   ,   b2k1T2 


2x, 


+ 


2x, 


A 


23,2 


b2ki       b2k{T2 

^34,2      =      ^24  H 1 ~ ^24,2 


Xd 


2x, 


3.  Third  order  approximation  in  y 

A  third  order  approximation  in  y  is  based  on, 


1 


(30) 


y(t-T)  =  y-Ty  +  l-T2y  -  -T^  ■ 

Similar  algebraic  steps  as  in  the  previous  two  approximations  result  in  the 
following  eigenvalue  problem, 


10       0  0       0 

0    10  0       0 

0    0    £33,3  0    £35,3 

0    0       0  10 

0    0    B53,3  0    #55,3  . 


v 

r 

y 

V2 
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0          0          10  0 

0          0          0          0  1 

^31,3      -432,3      -433,3      -434,3  -435,3 

110              0  0 


r 

y 


(31) 


-451,3     ^52,3      ^53,3      ^54,3      ^55,3  J    L  v2 

where,  v\  =  v,  v2  is  the  rate  of  change  of  v,  and  the  entries  of  the  generalized 
eigenvalue  problem  (31)  are  given  bellow.  Higher  order  approximations  in  T 
can  not  produce  any  usable  stability  results,  since  the  B  matrix  in  (31)  becomes 
singular. 


-43i,3    =    «n  +  Mi  - 


Xd 


-432,3     =     ai2  — 


bjkiT 

Xd 


-433,3      =      013  +  &1&2  + 


6ifciT; 
2xh 


-434,3      —      a14  + 


Xd 


-435,3 
^51,3 


b^T2 
2xd 

a2i  +  b2ki 


-452,3      =      G22  _ 


b2k1T 
xd 


^53,3      =      ^23  +  &2&2  + 


b2k1T 

Xd 


b2hT: 
2xd 


-454,3      =      &24  + 


Xd 


A 


55,3      — 


b2klT2 

2xd 
62fciT3 

6xd 

#35,3      =      £33,3 

b2hT3 


B 


33,3 


B 


53,3    — 


6xd 


+  1 


19 


B 


55,3      = 


b2k1Ti 
6xd 


4.  First  order  approximation  in  6 

Substitution  of  (26)  into  (28)  yields  the  following  matrix  system, 


1  0  0 

0  1  #23,4 

0  0  533,4     0 

0  0  0 


0  " 

l>] 

0 

V 

0 

r 

1  _ 

.  y  . 

0  0  10 

-421,4      ^22,4      -423,4      -424,4 

^31,4      -432,4      -433,4      ^34,4 

110  0 


where, 


121,4 


an  +  b\k\ 


&ifciT 

Xd 


^22,4      =      1x2  — 


Xd 


-423,4     =     013  +  b\k2  -  hkiT 


^24,4 


O14  + 


Xd 


-431,4      =      0-21  +  &2&1 


b2k1T 

Xd 


-432,4      =      0-22 


b2klT 

Xd 


V 

r 

y 


-433,4      =      ^23  +  b2k2  ~  b2k{T 

a  ,     &2&1 

-^34,4      =      Q-24  "I 

Xd 

#23,4      =      bYk2T 

#33,4     =     l  +  b2k2T 

5.  First  order  approximation  in  both  y  and  6 

Substitution  of  (27)  into  (28)  yields  the  following  matrix  system, 


10           0  0 

0      I?22,5  #23,5  0 

0      #32,5  #33,5  0 

0          0               0  1 


if; 
V 

r 

.  y 


0  0  10 

-421,5      -422,5      -423,5      -424,5 

-431,5      -432,5      -433,5      -434,5 

110  0 


V 

r 

y 


(32) 


,     (33) 
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where, 


„                      ,  ,  ,        Mi?i      hkiT2 
^21,5     =    an  +  dim 

bik\Ti       bikiT2 

Xd  Xd 

-423,5     =     ai3  +  0i«2  -  Ol«iTi  H 

6i/ci 


-422,5      =      Ol2  — 


-424,5      =      0-14  + 


Xd 

hk\Ti       b2kiT2 


-43i,5    =    0.21  +  b2ki  - 

Xd  xd 

b2k{Ti       b2k{T2 


-432,5      =      &22  — 


Xd  xd 

b2kiTlT2 


-433,5      =      «23  +  b2k2  ~  b2k{Ti  + 

Xd 

a                                 ,     ^l 
-4  34,5      =      ^24  H 

bikiTrT2 

#22,5      =      1 

#23,5      =      bik^ 

b2k1T1T2 

•D32,5      = 

Xd 

#33,5      =      1  +  b2k2Ti 


D.  EXACT  COMPUTATION  OF  TIME  LAG 

The  previous  analysis  using  Taylor  series  expansions  of  y(t  —  T)  breaks 
down  for  approximations  beyond  third  order,  as  the  corresponding  generalized 
eigenvalue  problem  becomes  singular.  Thus,  in  order  to  obtain  an  exact  com- 
putation of  the  stability  curves  and  check  the  validity  of  the  calculations,  a 
different  technique  will  be  performed  in  this  section.  This  technique  is  based 
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on  frequency  response  methods  and  it  utilizes  Nyquist's  criterion  for  stability 
[Friedland  (1986)]. We  write  the  system  of  equations  (11)  along  with  equations 
y  =  ip  +  v  and  8  —  k\ip  +  k2r  +  ~y(t  —  T),  in  the  Laplace  domain,  where 


Xd' 


.-Ts 


y(t-T)-^ye-'s,  (34) 

is  the  time  delay  operator. The  characteristic  equation  of  the  system  is, 

s4  +  Qi53  +  a2s2  +  a3s  +  q4  +  (/32s2  +  fas  +  /30)—e~Ts  =  0  ,         (35) 
where 

Oil      —      —0-12  -  «23  -  &2&2  , 

OL2      =      —^14  +  ^12^23  —  022^1^2  ~  ^22^13  +  a\2b2k2  ~  b2kj  _  «21  j 

0:3      =      —a24  —  Ol3a24  —  ^24^1^2  +  ^23014  +  01462^2  —  ^1^1^22  ~ 

^11^22  +  ^2^1^12  +  ^21^12  , 
#4      =      (221^14  +  ai2a24  —  ^22^14  _  ^11^24  +  b2k\a\4  —  bikid24  , 

02    =    ~hh  , 

01      =      -^13^2^1  +  ^1^1^23  -  &2^1  , 

Po      —     ^12^2^1  —  ^1^1^22  —  &2^1^H  +  bikiCL2l  , 

The  characteristic  equation  is  written  as, 

1  +  —  G(s)  =  0, 

Xd 


where 


s4  +  aiS"3  +  a2s*  +  Q3S  4-  Q4 
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is  the  open  loop  transfer  function  ,  and  ^-  denotes  the  effective  position  gain. 
For  stability,  we  can  utilize  the  Nyquist  criterion  which  states  that  the  critical 
value  of  Xd  can  be  computed  from, 

—  \(G(jcj)\  =  1        for         aigG(juj)  =  -n  (37) 

Xd 

where  j  denotes  the  imaginary  unit.  The  argument,  argC(jcj),  of(36)  is, 

nc    \             v,.      -i       01"            ,      -i     <*3^  ~  Qi^3  ,0Q, 

arg  G \j  (J)  =  -ujT  +  tan — --tan     — .  (38) 

The  set  of  equations  (37)  result  in  the  critical  visibility  distance, 


Xd 


$W  +  (A)  "  fou2)2 


(39) 


\  (asuj  —  aiuj3)2  +  (uj4  —  a2UJ2  +  a^)2 

The  value  of  u>  in(38)  such  that  argG(jco)  =  —  -k  is  the  value  of  the  phase 
crossover  frequency.  After  solving  for  the  phase  crossover  frequency,  the  critical 
value  of  Xd  is  obtained  by  direct  substitution  of  this  value  of  u  into  equation 
(39). 


E.  RESULTS  AND  DISCUSSION 

Typical  results  are  presented  in  Figures  3  through  6.  All  results  shown 
are  nondimensional  unless  otherwise  stated.  The  critical  value  of  Xd  versus  ujn 
for  zero  time  lag,  and  parametrized  for  different  values  of  the  damping  ratio 
£  is  shown  in  Figure  3.  Stability  is  ensured  for  values  of  Xd  greater  than  its 
critical  value.  It  can  be  seen  that  lower  values  of  ion  require  higher  values 
of  Xd  for  stability.    This  means  that  a  less  responsive  helmsman  will  need  a 
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Figure  3:  Critical  Xd  versus  wn  for  T<i  =  0  and  various  values  of  C 
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Figure  4:  Critical  x^  versus  u>n  for  C  =  0.8  and  various  values  of  T2 
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Figure  5:  Critical  x-d  versus  ujn  for  Q  —  0.8,  T<2  =  5sec  and  various  values  of  '1\ 
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Figure  6:  Critical  x<i  versus  un  for  C  =  0-8  and  zero  time  lag:  Channel  effects 


25 


longer  lookahead  distance  for  a  stable  operation.  Similar  conclusions  hold  for 
variations  in  the  damping  ratio  (.  In  this  case,  higher  values  of  £  correspond 
to  better  helmsman  response  which  requires  less  lookahead  distance. 

The  effect  of  time  lag  T2  is  shown  in  Figure  4.  Time  lags  are  in  seconds.  It 
can  be  seen  that  reasonable  amounts  of  time  lag  do  not  have  a  serious  effect  on 
stability,  at  least  in  a  linear  sense.  Of  course,  an  amount  of  time  lag  may  have 
a  serious  effect  on  the  transient  response  of  the  system  as  well  as  its  ability  to 
reject  an  external  disturbance.  This  can  only  be  established  by  a  systematic 
series  of  numerical  simulations.  Similar  conclusions  hold  for  non-zero  time 
lags  T\,  as  shown  in  Figure  5. 

Finally,  the  severe  destabilizing  effect  of  the  canal  is  demonstrated  by  the 
results  of  Figure  6  where  a  comparison  between  canal  and  open  water  results  is 
presented.  It  can  be  seen  that  an  order  of  magnitude  increase  in  the  lookahead 
distance  may  be  required  if  the  same  control  parameters  are  to  be  used  in  both 
open  water  and  a  canal. 
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IV.  NONLINEAR  ANALYSIS 


A.  INTRODUCTION 

It  can  be  numerically  confirmed  that  in  all  cases  of  stability  loss  of  the  pre- 
vious chapter,  one  pair  of  complex  conjugate  eigenvalues  of  the  corresponding 
eigenvalue  problem  crosses  transversally  the  imaginary  axis.  A  situation  like 
this  in  which  a  certain  parameter  is  varied  such  that  the  real  part  of  one  pair 
of  complex  conjugate  eigenvalues  of  the  linearized  system  matrix  crosses  zero, 
results  in  the  system  leaving  its  steady  state  in  an  oscillatory  manner.  This 
loss  of  stability  is  called  Hopf  bifurcation  and  generically  occurs  in  one  of  two 
ways,  supercritical  or  subcritical.  In  the  supercritical  case,  stable  limit  cycles 
are  generated  after  the  nominal  straight  line  motion  loses  its  stability.  The 
amplitudes  of  these  limit  cycles  are  continuously  increasing  as  the  parameter 
distance  from  its  critical  value  is  increased.  For  small  values  of  this  criticality 
distance  the  resulting  limit  cycle  is  of  small  amplitude  and  differs  little  from 
the  initial  nominal  state.  In  the  subcritical  case,  however,  stable  limit  cycles 
are  generated  before  the  nominal  state  loses  its  stability.  Therefore,  depend- 
ing on  the  initial  conditions  it  is  possible  to  diverge  away  from  the  nominal 
straight  line  path  and  converge  towards  a  limit  cycle  even  before  the  nominal 
motion  loses  its  stability.  This  means  that  in  the  subcritical  Hopf  bifurcation 
case  the  domain  of  attraction  of  the  nominal  state  is  decreasing  and  in  fact 
it  shrinks  to  zero  as  the  critical  point  is  approached.    Random  external  dis- 
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turbances  of  sufficient  magnitude  can  throw  the  vehicle  off  to  an  oscillatory 
steady  state  even  though  the  nominal  state  may  still  remain  stable.  After  the 
nominal  state  becomes  unstable,  a  discontinuous  increase  in  the  magnitude  of 
motions  is  observed  as  there  exist  no  simple  stable  nearby  attractors  for  the 
vehicle  to  converge  to.  Distinction  between  these  two  qualitatively  different 
types  of  bifurcation  is,  therefore,  essential  in  design.  The  computational  pro- 
cedure requires  higher  order  approximations  in  the  equations  of  motion  and  is 
the  subject  of  the  next  section. 


B.  DETAILED  CALCULATIONS 

The  nonlinear  equations  of  motion  are  written  as, 

ip    =    r  ,  (40) 

v     =    aiii>  +  ai2V  +  ai^r  +  auy  +  bi6' ,  (41) 

r    —    a2iip  +  a22v  +  a23r  +  a24y  +  b26'  ,  (42) 

y    =    sin  tp  +  v  cos  ip  ,  (43) 

where  the  control  law  is, 

i  y(t-T) 

8  =  kii)  +  k2r  +  fci  tan-1  — ,  (44) 

and,  including  saturation, 

£' =  £sattanh(-/-]   .  (45) 

\t>sat/ 

We  perform  a  Taylor  series  expansion  of  the  equations,  keeping  the  first  non- 
vanishing  nonlinear  terms.   Due  to  port/starboard  symmetry  in  the  problem 
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this  means  expansion  to  third  order  terms, 

simp  =ip  -  |V3  ,      cos?/;  =  1  -  \i\)2  ,  (46) 

'-°-ikf'  k  (47) 

8  =  k^  +  k2r  +  -iy(t  -  T)  -  -^y3{t  -  T)  ,  (48) 

where  for  consistency,  the  term  y(t  —  T)  is  approximated  by  its  first  order 
expansion  in  T, 

y(t-T)  =  y-Ty  =  y-TiP  +  \T^  -  Tv  +  \Tvi}j2  .  (49) 

Substitution  of  equations  (46)  to  (49)  into  equations  (40)  to  (45),  produces 
the  system, 

x  =  Ax  +  g(x),  (50) 

where  the  state  variables  vector  is, 

x  =  [ip,v,r,y}T  , 

A  is  the  linearized  matrix  at  equilibrium,  and  g(x)  contains  all  third  order 
terms. 

If  T  is  the  matrix  of  eigenvectors  of  A  evaluated  at  the  critical  point  Xd, 
the  linear  change  of  coordinates, 

x  =  Tz,      z  =  T_1x,  (51) 

transforms  system  (50)  into  its  normal  coordinate  form, 

z  =  T-1ATz4-T_1g(Tz)  .  (52) 
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At  the  Hopf  bifurcation  point,  matrix  T   *AT  takes  the  form, 


T_1AT  = 


0  -u>0  0  0 

w0       0  0  0 

0        0  p  0 

0        0  0  q 


where  luq  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 
remaining  two  eigenvalues  p  and  q  are  negative  (or  complex  conjugate  with 
negative  real  parts)  .  Therefore  in  the  new  system  of  coordinates  (51)  the 
dynamics  of  (50)  are  governed  by  a  reduced  two-dimensional  system  z\  and 
Z2  since  the  coordinates  23  and  24  correspond  to  the  eigenvalues  p  and  q  and 
are  asymptotically  stable.  For  values  of  Xd  close  to  the  bifurcation  point  Xd, 
matrix  T_1AT  is, 


T-iAT 


a'e  -{ujq  +  uj's)  0  0 

ojq  +  uj'e  ol'e  0  0 

0  0  p+p's         0 

0  0  0  q  +  q'e 


where,  s  denotes  the  criticality  difference, 


£  =  Xd  —  Xd   :♦• 

u  "critic 


(53) 


U)       = 


a     = 


p     = 


Q      = 


derivative  of  the  real  part  of  the  critical  eigenvalues  with  respect  to  e  , 

derivative  of  the  imaginary  part  of  the  critical  eigenvalues 

with  respect  to  e  , 

derivative  of  p  with  respect  to  e  , 

derivative  of  q  with  respect  to  e  . 


Due  to  continuity,  the  eigenvalues  p  +  p'e  and  q  +  q'e  remain  negative  for  small 
nonzero  values  of  e.  Therefore,  the  coordinates  23,  24  correspond  to  negative 
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eigenvalues  and  are  asymptotically  stable.  Center  manifold  theory  predicts 
that  the  relationship  between  the  critical  coordinates  z\,  z2  and  the  stable 
coordinates  23,  24  is  at  least  of  quadratic  order.  In  fact,  due  to  the  symmetry 
in  our  problem  the  relationship  is  cubic, 

21=e?(233,243),         22  =  C>(23,23). 

It  follows  that  because  of  this  order,  23,  24  do  not  influence  the  third  order 
Taylor  series  expansions,  and  can  be  dropped  from  the  equations.  Therefore, 
we  can  write  the  reduced  system  that  describes  the  essential  dynamics  of  (52) 
in  the  form, 

21    =    a'ezi  -  (u0  +  u'£)z2  + Fi(zi,z2)  ,  (54) 

i2     =     (ujq  +  ujr£)zi  +  a'sz2  +  F2(z1,z2)  ,  (55) 

where  F\,  F2  contain  the  third  order  terms, 

Fi(zi,z2)     =    ruzl  +  ri2ziz2  +  ri3ziz%  +  ruz%  ,  (56) 

^2(^1,22)     =    r2i23  +  r22z\z2  +  r2Zzizl  +  t24lz\  .  (57) 

The  coefficients  r#  are  computable  from  the  previous  Taylor  expansions,  and 
they  are  given  below. 

The  nonlinear  expansion  coefficients  r»j  that  are  utilized  in  the  definition 
of  the  cubic  stability  coefficient  K,  are  given  by  the  following: 

rn     =    n12£21  +  n13£3i  +  n14£41  +  n12dn  +  ™i3^2i  , 
T\2     =     n\2l22  +  7113^32  +  7114^42  +  n\2di2  +  rii3d22  , 
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ri3    =    ^12^23  +  "13^33  +  nu£43  +  ni2di3  +  ni3d23  , 

^14  =  ^12^24  +  ™13^34  +  ^14^44  +  "12^14  +  ni3d24  , 

1"21  =  n22hl  +  ™23^31  +  "24^41  +  "22^11  +  ™23^2X  , 

^22  =  ^22^22  +  ™23^32  +  ^24^42  +  ^22^12  +  ^23^22  , 

?~23  =  ^22^23  +  ™23^33  +  ^24^43  +  ^22^13  +  ^23^23  , 

T24  =  ^22^24  +  ^23^34  +  ^24^44  +  ^22^14  +  ^23^24  , 

where  we  denote, 

T=  [my]  ,  T-1  =  [nij]  ,  i,  j  =  1, . . .  ,4  . 

For  numerical  purposes  the  critical  eigenvector  of  T  must  be  normalized  so 
that  its  first  nonvanishing  coefficient  is  identically  equal  to  1.  The  coefficients 
£i:j  are  given  by, 

21  c2  c  2  o  2  c  2 

— —     =    <V^mnm2i  +  <Vwmnm2i  +  <-ty^rmnm3i  +  oT/,rrmiim31 
h 

+<5vV>ymiim4i  +  ^Vyymnm4i  +  ^wr"T-2im3i  +  ^^^21^31 

o  0  0  0 

+^W2/m21m41  +  6vyym2im4l  +  6rrym31m4i  +  5ryym3im41 

+^w"lll^21m31  +  <$i/wymnm21m41  +  ^^772117714177^3! 

+<5ury7n2i7n4im3i  +  S^^rrin  +  <5uwm2i  +  8rrrm2i  +  8yyymAl  , 

=      <V/^(mllm22  +  277T,nmi2m2l)  +  6vw(m12m21  +  2mxxmx2m22) 

+6^r  (mix  m32  +  2mnmi2^3i)  +  (Vrtm^mjn  +  2mnm3xm32) 
+<5^y(miim42  +  2mnmi2m4i)  +  ^viw(m4imi2  +  2miim4im42) 
+(^(77121777,32  +  2m3im2im22)  +  ^("122^3!  +  2m3im32m2i) 
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^22 
&1 


+SVVy(m2im42  4-  2m4im2im22)  +  Svyy(m22m41  +  2m4im42m2i) 

+6rry(m31m42  +  2m4im3im32)  +  ^7^(^132^41  +  2m4ira42ra3i) 

+6^vr[mum2im32  +  m3i(mnm22  +  mi2m2i)] 

+ S^vy  [mum  2im42  +  m4i(mnm22  +  m  1277121)] 

+S%pry[miim4imz2  +  ^31(^11^42  +  ^12^41)] 

+Svry[m2im,4ims2  +  7n3i(77i2i77i42  +  m22"T-4i)] 

+<5^V'tA(3m11mi2)  +  5wu(3m21m22)  +  <5rrr(3m31ra32)  +  5yyy(3m41m42)  , 

£ 

—  =    #vVf  (mi2m2i  +  2mnmi2m22)  +  fy*w(mnm22  +  2mi2m2iTn22) 
01 

+^#(^12^31  +  2mnmi2m32)  +  <Vr(mllm32  +  2mi2m3im32) 
+^y(mi2m4i  +  2mnmi2m42)  +  ^vyy  (m42mn  +  2777,12777,4177142) 

+  <5Wr(m22m31  +  2m21^22^32)  +  ^rr(m2im32  +  2m3im2277732) 

+<5wy(m227n4i  +  2m2im227n42)  +  <$7jyy(77i2i7n42  +  2777,41777,42777,22) 

+^TTj/(^32m41  +  2m42m3im32)  +  ^ryy(^31^42  +  2m4i77l42m32) 
+  ^ur.[m  1277222^31  +  77l32(777,i177l22  +  TTi^TT^l)] 
+^y[m12m227774i  +  m42(miim22  +  7711277121)] 

+<5^r.y[mi2m42m3i  +  m32  (777,11777,42  +  7711277141)] 
+6vry[m22,m42msi  +  m32{m2\m42  +  7712277141)] 
+8Tp^(3miiml2)  +  Svvv(3m2\m22)  +  6rrr(Smzim32)  +  6yyy(3m4im42)  , 

24  r2  c  2  c  2  c  2 

—  =    o-lp.l)}Vml2m22  +  o^vvmi2m22  +  d^rml2mz2  +  o^rr7ni27n32 
01 

+6-lp1pyml2m42  +  6^^7711277142  +  8vvrm22m32  +  <5wrm227n32 

o  9  9  9 

+<5uuy7n2277i42  +  bVyyrn22'm42  +  Srrym32m42  +  6ryym32m42 
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+6^vrmi2m22'ms2  +  <Vvym  12^22^42  +  5'0rym12m42'7l32  +  ^ury^22^42^32 

4i         4i 

62  61 

^32  ^22 

62  &1 

^33  ^23 

b2  h 

h4  ^24 

b2  h 

1     2    /  1 

^41  =    -2mn  (m2i  + 3mn 

£42  =     -mn  ( mi2m2i +  -mnm22  + -mi2mn  J   , 

^43  =      -mi2  f  mnm22  +  -7^1277121  +  -mi2mii  J     , 

£44  =      --77li2  (m22  +  -77112 J     • 


The  coefficients  da  are  given  by, 


dn     =    —[cu(Iz  -  Nr)  +  c2i{Yr  -  mxG)]  , 

i-J  vr 
d\2      =      —  [Cl2(lz  ~  Nr)  +  C22{Yr  ~  mXG)]  , 

J-JyT 

dl3      =      —[Ci3(lz  ~  Nf)  +  C23(Yf  -  mxG)]  , 
du     =     — —  [cu(Iz  -  Nr)  +  c24(Yf  ~  mxG)}  , 

L>vr 

d2\     =    -r—[cu{Ny  -  mxG)  +  c2i(m  -  Y*)]  , 

L)VT 

d22    =     r=—[ci2(Ni  -  mxG)  +  c22(m  -  ¥&)]  , 
u  vr 

^23     =     jr-[ciz{Ny  -  mxG)  +  C2z{m  -  Yv)}  , 
d2A    =    -fr-[cu(Ni  -  mxG)  +  c24(m  -  Y*)]  , 
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where 


o  J.  o  J.  o 

en    =    -Yvvvm21  +  7:Yvrr™>3lm,2i  4  -YTVVm2lmz\ 

1  3  1  3 

+     -}'WV'mll  +  7^yyj/m4l  > 

c12    =     -Yvvv3mllm22  4  -zYvrr{m2zlm22  +  2m3im32m2i) 
4     7Yrvv{™>2imz2  4  2m3im2im22)  4-  -Y^^m^m^ 

1  2 

+      T^yyy3m41m42  , 

C13    =     -Vt;tw3m22m2i  4  -ywr.(m32m2i  4  2m31m32m22) 

D  Z 

+      -^rw(^22m31  +  2m32m2im22)  4-  ->V'V>V'3rn12rnll 
1  2 

4  ^YyyySm  ^2'm  41     1 

X  O  J.  Q  JL  O 

C14      =      -yrwvm22  4  -^rr^32m22  +  -Yrvvm22m32 
O  Z  2 

4     ~Y^^ml2  4  ^YyyyrnA2  > 

o  J.  o  J.  o 

C21     =     -Nvvvm2l  +  -iVwrm31m2i  4  -Nrvvm21m3i 

0  z  z 

1  3  1  3 

4     —N^^rrin  +  -Nyyym4l  , 
C22    =     -Ntw«3m2im22  4  -^^(^3^22  4  2m31m32m2i) 

D  Z 

4     -ArrVV(^2im32  4  2m31m2im22)  4  -A/'^v3miimi2 

Z  D 

1  2 

4     ~A/'yyy3m41m42  , 
C23    =     -NvvvZm22'm,2\  4  -N^m^r^i  4  2m3imz2m22) 

0  Z 


i/V_....frr2   -  '    «- 

2 

1- 


4     -Nrvv{m\2mz\  4  2m32m2im22)  4  -N^^m^mn 


4     -Nyyydm^m^  , 

3  2  2 

C24     =     -Nvvvm22  +  -Nvrrm32m22  +  -Nrvvm22m32 

1  3        1  3 

4     -N^^ml2  4  -Nyyym42  . 

0  D 
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The  coefficients  in  a  third  order  expansion  of  the  control  law  are  defined 

by 

A;iT3 


1                               T       k\ 
S^ijjv    =    —  -To- {k[)  k'2  +  0.5/ci 1 


''sat  *d  ■*•' 


rf 


dxpvv      —  co        IV    2/       '  3 


s 


sat 

0Sat 

X. Lt'1-2 

uyirr      —  cO        12  ' 


82 


sat 


/\2/tl 


fei     r2^ 


°sat  ^d  xd 

1       ,  fc?         T/Ci 


-iri 

°sat       x  - 


<->i/>yy      —  j-2        1  ^2 

'S£ 

Cr      =      ~T2     v    2/    ^2  j 
^sat 

X Lt-'t-2 

uurr      —  <-2     A'2"'2  ' 

dsat 


1    Y£/\2fcl  ^l^2 

sat 


<$wy       —  .2     (*2)     _ 

°sat  xd  xd 

°vyy      —  r2    ^2     9  ^     _3     ' 

°sat       xd  xd 

1       ofci 

°rry      —  <-2        2  ' 

^sat        ^-d 

1        /c2 

°ryy      —  r2    *-2     2   ' 

°sat       ^-d 

Vxjjvr      ==  T2     ^Ki'^2'^2  > 

<->sat 

1  ,     ,*i  T2/d 

^Vwy     —     —  T2    ^"-1*2  " 


^sat  ^d  Xc 

Oi/jry      ==      —  T2     ^"-l"-2  j 

<->sat  ^d 

1  ,       fei 

<5s2at  *d 
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,  1    1 ,  ,,s3  ,   W      fc,T3 


ysat  "-"  uxd  oxd 


1         !,_..*  fclT 


3 


0--+  O  OX  j 


'sat 
1     1 


X  I_  .2  1.3 

yrrr        —  c2       Q      2    ' 


sat 

1    \k\       2/ci 


dent  O  Ij  OIj 


and 
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/C-]         —         /C^  rCj  - 

%d 

T 
fc2    =     —  k\       . 

Xd 


C.  NONLINEAR  COEFFICIENT  K 

If  we  introduce  polar  coordinates  in  the  form, 

z\  =  RcosO  ,      22  =  -Rsinfl  ,  (58) 

we  can  use  (54)  and  (55)  to  produce  an  equation  describing  the  rate  of  change 
of  the  radial  coordinate  R, 

R  =  a'eR  +  V{6)Rz  ,  (59) 

and  a  similar  equation  in  the  rate  of  change  of  the  angular  coordinate  0, 

6  =  uj0  +  uj'e  +  Q{9)R2  .  (60) 

The  system  of  equations  (59)  and  (60)  contains  two  variables,  one  of  which, 
i?,  is  slowly  varying  in  time,  whereas  the  other  one,  6  is  a  fast  variable.  Then, 
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equation  (59)  can  be  averaged  over  one  cycle  in  0  to  produce  an  equation  with 
constant  coefficients  and  similar  stability  properties, 

R  =  a'eR  +  K,RZ  ,  (61) 

where, 

1     f2n 
K  =  —         V{9)  d9  =  |  (3rn  +  na  +  r22  +  3r24)  .  (62) 

Z7T  70 

Nontrivial  equilibrium  solutions  of  (61)  correspond  to  limit  cycles  in  the 
original  system.  The  nontrivial  equilibrium  of  (61),  Rq,  is  given  by, 


a's 


Ro  =  \J-^-  (63) 

In  our  problem,  since  the  trivial  equilibrium  gains  its  stability  for  xj  >  ^dcritical, 
the  coefficient  a'  is  always  negative.  Therefore,  it  can  be  seen  from  (63)  that 
a  limit  cycle  will  exist  provided  K  and  e  have  the  same  sign.  The  stability 
properties  of  this  limit  cycle  can  be  determined  by  linearization  of  (61)  around 
Rq.  The  linearized  system  is, 

R  =  -2a'e(R  -  R0)  ,  (64) 

and  its  eigenvalue  is, 

(3  =  -2a  e  .  (65) 

We  can  see  that  the  Floquet  exponent  (65)  is  negative  if  e  is  negative,  which 
means  that  K.  must  be  negative.  Therefore,  location  and  stability  of  limit 
cycles  depends  entirely  on  the  cubic  coefficient  JC  which  is  computable  from 
(62).  We  can  summarize  our  findings  as, 
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Figure  7:  Nonlinear  coefficient  K,  versus  un  for  <5sat  =  0.4  and  various  values  of  C, 

•  If  K  <  0  then  limit  cycles  exist  for  e  <  0  (xj  <  £<*„.„.;„,)  and  they  are 
stable. 

•  If  K  >  0  then  limit  cycles  exist  for  e  >  0  (rc^  >  £dcritical)  and  they  are 
unstable. 


D.  RESULTS  AND  DISCUSSION 

Typical  results  in  terms  of  the  nonlinear  stability  coefficient  /C  are  shown  in 
Figures  7  through  10.  The  general  conclusion  from  this  series  of  graphs  is  that 
the  bifurcations  are  all  strongly  subcritical.  This  means  that  any  linearized 
stability  results  should  be  viewed  with  extreme  caution.    Limit  cycles  exist 
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Figure  8:  Nonlinear  coefficient  K,  versus  uin  for  C  =  0.8  and  various  values  of  <5sat 
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Figure  9:    Nonlinear  coefficient  K.  versus  un  for  Q  =  0.8,  <5sat  =  0.4,  T\  —  0,  and 
various  values  of  T2  (sec) 
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Figure  10:  Nonlinear  coefficient  /C  versus  wn  with  and  without  channel  effects 

before  stability  in  the  linear  sense  is  lost  and  a  self  sustained  oscillation  in  the 
system  may  develop  as  a  result  of  an  external  disturbance  even  if  the  nominal 
equilibrium  state  is  still  stable.  The  bifurcations  become  more  subcritical; 
i.e.,  /C  is  more  positive,  for  smaller  values  of  £,  as  Figure  7  shows.  Figure  8 
demonstrates  the  effect  that  the  saturation  limit  8sat  has  on  the  value  of  K. 
Higher  values  of  53at,  although  are  not  related  to  the  critical  value  of  Xd  result 
in  significant  changes  in  the  value  of  K,.  The  general  trend  is  that  higher  <5sat 
results  in  less  subcritical  bifurcations  as  evidenced  by  the  overall  decrease  of 
/C.  Figure  9  shows  the  effect  of  non-zero  time  lag  on  the  nonlinear  stability 
coefficient.  It  can  be  seen  that,  like  the  linear  stability  results,  the  effect  is 
minimal.    Finally,  Figure  10  shows  the  canal  effect  on  K.    This  figure  was 
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produced  for  zero  time  lag,  (  =  0.8,  and  <5sat  =  0.4.  It  can  be  seen  that  the 
existence  of  a  canal  causes  the  bifurcations  to  be  much  more  subcritical  than 
the  open  water  case.  This  demonstrates  the  severe  destabilizing  effect  that  the 
canal  introduces  in  both  the  linearized  and  the  nonlinear  analysis. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  comprehensive  study  of  linear  and  nonlinear  stability 
properties  of  straightline  motion  of  surface  ships  in  confined  waters.  The 
classical  maneuvering  equations  of  motion  incorporating  canal  suction  effects, 
were  coupled  with  appropriate  navigation,  gauidance,  and  control  laws  in  order 
to  mimic  the  helmsman's  behavior.  The  main  conclusions  from  this  study  can 
be  summarized  as  follows: 

1.  There  exists  a  critical  preview  distance  for  straightline  positional  stabil- 
ity. For  values  less  than  the  critical  distance,  the  system  is  unstable. 

2.  Including  canal  effects,  the  critical  preview  distance  may  be  an  order  of 
magnitude  higher  than  in  open  water.  If  the  same  preview  distance  is  to 
be  used  in  both  cases,  the  corresponding  control  law  for  canal  maneuver- 
ing must  be  considerably  more  responsive  than  in  open  waters. 

3.  The  critical  preview  distance  is  monotonically  decreasing  for  increasing 
control  law  responsiveness.  This  means  that  in  order  to  accomodate 
smaller  values  of  the  preview  distance,  more  responsive  control  laws  are 
required  . 

4.  Physically  realizable  time  lags  do  not  seem  to  have  a  significant  effect  on 
the  values  of  the  preview  distance. 

5.  As  the  preview  distance  becomes  less  than  its  critical  value,  one  pair  of 
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complex  conjugate  eigenvalues  of  the  linearized  system  matrix  crosses 
the  imaginary  axis.  This  corresponds  to  a  bifurcation  to  periodic  so- 
lutions (Hopf  bifurcation)  and  the  system  exhibits  oscillatory  behavior, 
also  known  as  limit  cycles. 

6.  Higher  order  approximations  in  the  equations  of  motion  were  utilized  in 
order  to  assess  the  stability  of  the  resulting  limit  cycles.  It  was  found 
that  in  all  cases,  the  limit  cycles  were  unstable.  This  has  the  following 
implications: 

(a)  It  is  possible  for  the  system  to  lose  its  stability  even  before  the  critical 
preview  distance  is  crossed.  This  means  that  all  linearized  stability 
analysis  results  should  be  viewed  with  extreme  caution. 

(b)  As  the  critical  preview  distance  is  crossed,  it  is  expected  that  the 
system  will  develop  limit  cycles  of  large  amplitude.  This  clearly 
presents  a  dangerous  situation  which  should  be  avoided  in  practice, 
by  appropriate  changes  in  the  design  parameters. 

Recommendations  for  further  research  include  the  following: 

1.  Simulation  studies  in  order  to  verify  the  limit  cycle  behavior. 

2.  Study  of  different  ship  characteristics  and  canal  geometry. 
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APPENDIX 


The  following  is  a  list  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN.  Complete  printouts  of  the 
programs  follow.  Standard  eigenvalue/eigenvector  numerical  analysis  subrou- 
tines are  required  for  all  programs. 

•  THESIS1.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lag  Ti- 

•  THESIS2.FOR 

Calculation  of  critical  Xd-  Second  order  approximation  for  time  lag  T2. 

•  THESIS3.FOR 

Calculation  of  critical  Xd-  Third  order  approximation  for  time  lag  T2. 

•  THESIS4.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lag  T\. 

•  THESIS5.FOR 

Calculation  of  critical  Xd-  First  order  approximation  for  time  lags  T\  and 
T2. 

•  HOPF.FOR 

Calculation  of  the  nonlinear  cubic  stability  coefficient  /C.    Requires  the 
output  of  one  of  the  previous  programs  as  its  input. 
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C     PROGRAM  THESIS 1. FOR  (Time  Delay-lst  Order  Approx.  T2) 

C     BIFURCATION  ANALYSIS 

C 

PROGRAM  THESIS 1 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,L,NR,NV,NDELTA,NPSI ,NY,IZ,MASS, 
&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,FV1(4) , IV1(4) ,ZZZ(4,4) ,WR(4) ,WI(4) 
C 

OPEN  (11,FILE='BIF1.RES') 

OPEN  (12,FILE='BIF2.RES') 

OPEN  (13,FILE='BIF3.RES') 

open  (15, f ile='eig.resl' ) 


c 

c 

Vehicle 

;  Parameters: 

IZ 

=0.0 

L 

=528 

RHO 

=  1.94 

XG 

=0.0 

MASS 

=0.0088 

U 

=  1.0 

c 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0.00456 

YV 

=-0.01434 

YPSI 

=  0.01400 

YY 

=  0.02000 

YDELTA=  0.00278 

NRDOT 

=-0.00115 

NVDOT 

=  0.00000 

NR 

=-0.00296 

NV 

=-0.00460 

NPSI 

=  0.01000 

NY 

=-0.00250 

NDELTA=-0. 00139 

WRITE 

(*,1001) 

READ 

(*,*)    WNMIN.WNMAX.IWN 

WRITE 

(*,1002) 

READ 

(*,*)     XDMIN,XDMAX,IXD 

WRITE 

(*,1003) 

READ 

(*,*)     ZETA 

WRITE 1 

>,1100) 

READ  ( 

[*,*)    TL 

TL=TL*U/L 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
&     (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 
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AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( (IZ-NRDOT) * (YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) * (NR-MASS*XG)) /DVR 
AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG)*NY) /DVR 
AA24=( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 
BB1  =( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 
BN0M=BB2 
CN0M=AA21 
EPS  -l.D-5 
ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 
WN=WNMIN+(I-1)*(WNMAX-WNMIN)/(IWN-1) 
Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 
K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+(J-1)*(XDMAX-XDMIN)/(IXD-1) 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
&    AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
write (15 ,*)DEOS,XD)WN 

IF  (J.GT.l)  GO  TO  10 
DEOSOO=DEOS 
XDOO  =XD 
LL=0 
GO  TO  2 
10     DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
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DEOSN=DEOSNN 
6     XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 

DEOSR=DEOSN 

XD=(XDL+XDR)/2.D0 
C 

CALL  LINEAR(K1)K2,XD,AA11,AA12,AA13,AA14,AA21, 
&     AA22,AA23,AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DE0SM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS(XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 
5     IF  (PRR.GT.O.DO)  STOP  3200 

XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 
4     LLL=10+LL 

WRITE  (LLL,*)  XD.WN 
3     XDOO=XDNN 

DEOSOO=DEOSNN 
2   CONTINUE 
1  CONTINUE 
C 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless)  ' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL  (dimensional)') 
2001  FORMAT  (215) 


48 


END 
C 

SUBROUTINE  DSTABL (DEOS , WR , WI , OMEGA) 
C 

C     EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR(4),WI(4) 
DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 
DEOS=WR(I) 
IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 
C 

SUBROUTINE  LINEAR(K1,K2,XD,AA11 ,AA12,AA13,AA14,AA21, 
&  AA22,AA23,AA24,BB1,BB2,A,TL) 
C 

C     FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation 
C 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  K1.K2 

DIMENSION  A (4, 4) 

A(l,l)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(l,4)=O.ODO 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 

A (2 , 2) =AA12-BB1*K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A(2,4)=AA14+BB1*K1/XD 

A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 

A(3 , 2) =AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A(3,4)=AA24+BB2*K1/XD 

A(4,l)=1.0DO 

A(4,2)=1.0DO 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 
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C     PROGRAM  THESIS2.F0R  (Time  Delay-2nd  Order  Approx  T2) 

C     BIFURCATION  ANALYSIS 

C 

PROGRAM  THESIS2 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 
&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,WR(4) ,WI(4) 


OPEN  (11,FILE='BIF1.RES') 

OPEN  (12,FILE='BIF2.RES') 

OPEN  (13,FILE='BIF3.RES') 

open  (15,f ile='eig.resl') 


c 

c 

Vehicle  Paramete: 

IZ 

=0.0 

L 

=528 

RHO 

=1.94 

XG 

=0.0 

MASS 

=0 . 0088 

U 

=  1.0 

c 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0.00456 

YV 

=-0.01434 

YPSI  =  0.01400 
YY  =  0.02000 
YDELTA=  0.00278 

NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV    =-0.00460 

NPSI  =  0.01000 
NY    =-0.00250 

NDELTA=-0.00139 
WRITE  (*,1001) 

READ  (*,*)    WNMIN,WNMAX,IWN 

WRITE  (*,1002) 

READ  (*,*)     XDMIN,XDMAX,IXD 

WRITE  0,1003) 
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READ   (*,*)    ZETA 
WRITE (*, 1100) 
READ  (*,*)  TL 
TL=TL*U/L 

DVR  = (IZ-NRDOT)* (MASS-YVDOT) - 
&     (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11= ( ( IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *  YV) /DVR 
AA13= ( (IZ-NRDOT) * (YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) *(NR-MASS*XG) ) /DVR 
AA14= ( (IZ-NRDOT) *YY+(YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB1  = ( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 

DO  1  I=1,IWN 

WRITE  (*,2001)  I.IWN 
WN=WNMIN+(I-1)*(WNMAX-WNMIN)/(IWN-1) 
Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 

DO  2  J=1,IXD 

XD=XDMIN+(J-1)*(XDMAX-XDMIN)/(IXD-1) 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
&    AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERH) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
write ( 15 ,*)DEOS,XD,WN 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO   =XD 
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LL=0 
GO  TO  2 
10     DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6     XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR(K1,K2,XD,AA11)AA12,AA13,AA14,AA21, 
&     AA22,AA23,AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
DEOSM=DEOS 
XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 
IF  (PRL.GT.O.DO)  GO  TO  5 
XDO=XDL 
XDN=XDM 
DEOSO=DEOSL 
DEOSN=DEOSM 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDL-XDM) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
GO  TO  4 
5     IF  (PRR.GT.O.DO)  STOP  3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
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4     LLL=10+LL 

WRITE  (LLL,*)  XD.WN 
3     XDOO=XDNN 

DE0S0O=DEOSNN 
2   CONTINUE 
1  CONTINUE 
C 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless) ' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless) ' ) 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL  (dimensional)') 
2001  FORMAT  (215) 
END 
C 

SUBROUTINE  DSTABL (DEOS , WR , tf I , OMEGA) 
C 

C     EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR(4),WI(4) 
DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 
DEOS=WR(I) 
IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 
C 

SUBROUTINE  LINEAR(K1 ,K2,XD,AA11 ,AA12,AA13,AA14,AA21, 
&        AA22,AA23,AA24,BB1,BB2,A,TL) 
C 

C     FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation) 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DOUBLE  PRECISION  K1.K2 
DIMENSION  A (4, 4) 
A(1,1)=0.0D0 
A(1,2)=0.0D0 
A(1,3)=1.0D0 
A(1,4)=0.0D0 

A(2 , 1) = (AA11*XD+BB1*K1*XD-BB1*K1*TL) / (XD-0 . 50D0*BB1*K1*TL*TL) 
A(2,2)=(AA12*XD-BB1*K1*TL)/(XD-0.50D0*BB1*K1*TL*TL) 
A(2,3)=(AA13*XD+BB1*K2*XD+0.50D0*BB1*K1*TL*TL)/ 
&        (XD-0.50D0*BB1*K1*TL*TL) 
A(2,4)=(AA14*XD+BB1*K1)/(XD-0.50D0*BB1*K1*TL*TL) 
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A (3 , 1) =AA21+BB2*K1-BB2*K1*TL/XD+ (BB2*K1*TL*TL/ (2 . ODO*XD) ) *A(2 , 1) 
A(3,2)=AA22-BB2*K1*TL/XD+(BB2*K1*TL*TL/(2.0D0*XD))*A(2,2) 
A (3 , 3) =AA23+BB2*K2+ (BB2*K1 *TL*TL/ (2 . ODO*XD) ) + 
k  (BB2*K1*TL*TL/(2.0*XD))*  A (2, 3) 

A (3 , 4) =AA24+BB2*K1/XD+ (BB2*K1*TL*TL/ (2 . ODO*XD) ) *A(2 , 4) 
A(4,l)=i.0D0 
A(4,2)=1.0D0 
A(4,3)=0.0D0 
A(4,4)=0.0D0 
RETURN 
END 
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C     PROGRAM  THESIS3.F0R  (Time  Delay-3rd  Order  Approx  T2) 

C     BIFURCATION  ANALYSIS 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI ,NY, 
&  NRDOT , NVDOT 

DIMENSION  A(5,5) ,B(5,5) ,FV1(5) ,IV1(5) ,ZZZ(5,5) ,ALFR(5) , 
&  ALFI(5),BETA(5),WR(5),WI(5) 

C 

OPEN  (11,FILE='BIF1.RES') 

OPEN  (12,FILE='BIF2.RES') 

OPEN  (13,FILE='BIF3.RES') 


c 

c 

Vehicle  Parame 

IZ 

=0.0 

L 

=528 

RHO 

=  1.94 

XG 

=0.0 

MASS 

=0.0088 

U 

=  1.0 

YRDOT  =  0.00000 
YVDOT  =-0.00912 
YR  =+0.00456 
YV  =-0.01434 
YPSI  =  0.01400 
YY  =  0.02000 
YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV  =-0.00460 
NPSI  =  0.01000 
NY  =-0.00250 
NDELTA=-0.00139 


WRITE  (*,1001) 

READ   (*,*) 

WNMIN,WNMAX,IWN 

WRITE  (*,1002) 

READ   (*,*) 

XDMIN,XDMAX,IXD 

WRITE  0,1003) 

READ   (*,*) 

ZETA 

WRITEO.1100) 

READ  (*,*)  TL 

TL=TL*U/L 

DVR  =(IZ-NRDOT)*(MASS-YVDOT)- 
l  (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 
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AM  1= ( ( IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( (IZ-NRDOT) * ( YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) * (NR-MASS*XG) ) /DVR 
AA14= ( (IZ-NRDOT) *YY+ (YRDOT-MASS+XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB1  = ( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  -l.D-5 

ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 
WN=WNMIN+(I-1)*(WNMAX-WNMIN)/(IWN-1) 
Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 
K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+( J-l) * (XDMAX-XDMIN) / (IXD-1) 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 
&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (5 , 5 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , IERR) 
DO  11  IJE-1,5 

WR(IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
11     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

IF  (J.GT.l)  GO  TO  10 
DEOSOO=DEOS 
XDOO  =XD 
LL=0 
GO  TO  2 
10     DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 
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IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DE0S0=DE0S00 
DEOSN=DEOSNN 
6     XDL=XDO 
XDR=XDN 
DE0SL=DE0S0 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
&        AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (5 , 5 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , IERR) 
DO  12  IJE=1,5 

WR ( IJE) =ALFR ( IJE) /BETA ( IJE) 
WI ( IJE) =ALFI( IJE) /BETA (IJE) 
12     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS (XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 

GO  TO  4 
5     IF  (PRR.GT.O.DO)  STOP  3200 

XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS(XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 
4     LLL=10+LL 

WRITE  (LLL,*)  XD.WN 
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3     XDOO=XDNN 

DEOSOODEOSNN 
2   CONTINUE 
1  CONTINUE 


1001  FORMAT  ( 

1002  FORMAT  ( 

1003  FORMAT  ( 
1100  FORMAT  (: 
2001  FORMAT  (215) 

END 


ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless) ' ) 
ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 
ENTER  DAMPING  RATIO') 
ENTER  TIME  LAG  TL  (dimensional)') 


SUBROUTINE  DSTABL (DEOS , WR , WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(5),WI(5) 
DE0S=-1.0D+20 
DO  1  1=1,5 

IF  (WR(I).LT.DEOS)  GO  TO  1 

DEOS=WR(I) 

IJ=I 
L  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 

SUBROUTINE  LINEAR(K1 ,K2,XD,AA11 ,AA12,AA13,AA14,AA21,AA22,AA23, AA24 
&  ,BB1,BB2,A,B,TL) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2 

DIMENSION  A(5,5),B(5,5) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(1,5)=0.0D0 

A(2,1)=0.0D0 

A(2,2)=0.0D0 

A(2,3)=0.0D0 

A(2,4)=0.0D0 

A(2,5)=1.0D0 

A(3,1)=AA11-(BB1*K1*TL/XD)+BB1*K1 

A(3,2)=AA12-BB1*K1*TL/XD 

A (3 , 3) =AA13+BB1*K2+ (BB1*K1*TL*TL/ (2 . DO*XD) ) 

A(3,4)=AA14+BB1*K1/XD 

A(3,5)=-1.D0+(BB1*K1*TL*TL/(2.D0*XD)) 

A(4,1)=1.0D0 
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A 

A 

,2 

)=1.0D0 

A 

A 

,3 

)=0.0D0 

A 

A 

,4 

)=0.0D0 

A 

A 

,5 

>=0.0D0 

A 

[s 

,1 

>=AA21+BB2*K1-BB2*K1* 

fL/XD 

A 

:5 

2 

)=AA22-BB2*K1*TL/XD 

A 

:5 

,3 

) =AA23+BB2*K2+ (BB2*K1*TL*TL/ (2 . DO*XD)  ) 

A 

[5 

4 

)=AA24+BB2*K1/XD 

A 

:5 

5 

>=(BB2*K1*TL*TL/(2.D0 

*XD)) 

B 

;i 

i: 

)=1.0D0 

B< 

;i 

2] 

1=0. ODO 

B< 

:i 

3: 

I =0 . ODO 

B( 

;i 

4: 

1=0. ODO 

B( 

:i 

5: 

1=0. ODO 

B( 

2 

i: 

1=0. ODO 

B( 

2 

2] 

1=1. ODO 

B( 

[2 

3] 

1=0. ODO 

B( 

2 

4: 

»=0.0D0 

B( 

2 

5: 

1=0. ODO 

B( 

3 

i: 

=0 . ODO 

B( 

3 

2] 

=0 . ODO 

B( 

3 

3: 

=(BBl*Kl*TL*TL*TL/(6 

.DO*XD)) 

B( 

3 

4; 

=0 . ODO 

B( 

3 

5; 

=(BBl*Kl*TL*TL*TL/(6 

.DO*XD)) 

B( 

4 

1: 

=0.0D0 

B( 

4 

2; 

=0 . ODO 

B( 

4 

3; 

=0 . ODO 

B( 

4 

4; 

=1.0D0 

B( 

4 

5; 

=0 . ODO 

B( 

5 

1: 

=0 . ODO 

B( 

5 

2; 

=0 . ODO 

B( 

5 

3; 

=1.0D0+(BB2*K1*TL*TL> 

*TL/(6.D0*XD)) 

B( 

5 

4; 

=0 . ODO 

B( 

[5 

5: 

= (BB2*K1*TL*TL*TL/ (6 

.DO*XD)) 

RI 

:turi 

I 

El 

fD 
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C     PROGRAM  THESIS4.F0R  (Time  Delay-lst  Order  Approx  in  delta  ,T1) 

C     BIFURCATION  ANALYSIS 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,L,NR,NV,IZ,MASS ,NDELTA,NPSI ,NY, 
&  NRDOT.NVDOT 

DIMENSION  A(4,4),B(4,4) ,FV1(4) , IV1(4) ,ZZZ(4,4) , ALFR(4) , 
&  ALFI(4),BETA(4),WR(4),WI(4) 

C 

OPEN  (11,FILE=,BIF1.RES') 

OPEN  (12,FILE='BIF2.RES') 

OPEN  (13,FILE='BIF3.RES') 
C 
C    Vehicle  Parameters: 

IZ    =0.0 

L     =528 

RHO   =1.94 
C     G     =32.2 

XG    =0.0 

MASS  =0.0088 

U     =3.0 


YRDOT  =  0.00000 
YVDOT  =-0.00912 
YR  =+0.00456 
YV    =-0.01434 

YPSI  =  0.01400 
YY    =  0.02000 

YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV    =-0.00460 

NPSI  =  0.01000 
NY    =  -0.00250 

NDELTA=-0. 00139 


WRITE 

(* 

1001) 

READ 

(* 

*) 

WRITE 

(* 

1002) 

READ 

(* 

*) 

WRITE 

(* 

1003) 

WNMIN,WNMAX,IWN 
XDMIN.XDMAX.IXD 
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READ   (*,*)    ZETA 
WRITE(*,1100) 
READ  (*,*)  TL 
TL=TL*U/L 

DVR  = (IZ-NRDOT)* (MASS-YVDOT) - 
&     (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11= ( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI ) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) * (YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) * (NR-MASS*XG)) /DVR 
AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG) * NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB1  =( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2  =( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 
BN0M=BB2 

EPS  =l.D-5 
ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 

WN=WNMIN+(I-1)*(WNMAX-WNMIN)/(IWN-1) 

K1=-WN**2.D0/BN0M 

K2=- (AN0M+2 . DO*ZETA*WN) /BNOM 

DO  2  J=1,IXD 

XD=XDMIN+(J-1)*(XDMAX-XDMIN)/(IXD-1) 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 
&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , IERR) 
DO  11  IJB-1,4 

WR ( IJE) =ALFR ( I JE) /BETA ( I JE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
11     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

IF  (J.GT.l)  GO  TO  10 

DEOSOO=DEOS 

XDOO  =XD 
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LL=0 
GO  TO  2 
10     DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6     XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
&        AA23,AA24,BB1,BB2,A,B,TL) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , IERR) 
DO  12  IJE-1,4 

WR(IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
12     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 
XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 
IF  (PRL.GT.O.DO)  GO  TO  5 
XDO=XDL 
XDN=XDM 
DEOSO=DEOSL 
DEOSN=DEOSM 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDL-XDM) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
GO  TO  4 
5     IF  (PRR.GT.O.DO)  STOP  3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
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IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS(XDM-XDR) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
LLL=10+LL 

WRITE  CLLL,*)  XD.WN 
XD00=XDNN 
DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 


1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless) ' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL  (dimensional)') 

2001  FORMAT  (215) 
END 

SUBROUTINE  DSTABL(DEOS,WR,WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 
DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 

DEOS=WR(I) 

IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS( OMEGA) 
RETURN 
END 

SUBROUTINE  LINEAR (K1,K2,XD,AA 11 ,AA12,AA13,AA14,AA21, AA22,AA23,AA24 
&  ,BB1,BB2,A,B,TL) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2 

DIMENSION  A(4,4) ,B(4,4) 

A(1,1)=O.ODO 

A(l,2)=0.OD0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 
A(2 , 2) =AA12-BB1*K1*TL/XD 
A(2,3)=AA13+BB1*K2-BB1*K1*TL 
A(2 , 4) =AA14+BB1*K1/XD 
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3,1 
3,2 
3,3 
3,4 

4,1 
4,2 
4,3 
4,4 

1,1 
1,2 
1,3 
1,4 

2,1 
2,2 
2,3 
2,4 

3,1 
3,2 
3,3 
3,4 

4,1 
4,2 
4,3 
4,4 


=AA21+BB2*K1-BB2*K1*TL/XD 
=AA22-BB2*K1*TL/XD 
=AA23+BB2*K2-BB2*K1*TL 
=AA24+BB2*K1/XD 

=1.0D0 
=1.0D0 
=0 . 0D0 
=0 . ODO 

=1.0D0 
=0.0D0 
=0.0D0 
=0 . ODO 

=0 . ODO 
=1.0D0 
=BB1*K2*TL 
=0 . ODO 

=0 . ODO 
=0 . ODO 

=1.0D0+(BB2*K2*TL) 
=0 . ODO 

=0.0D0 
=0 . ODO 
=0.0D0 
=1.0D0 


RETURN 
END 
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C     PROGRAM  THESIS5.F0R  (Time  Delay-lst  Order  Approx  in  delta  &  y  ie.  in  Tl  &  T2)) 

C     BIFURCATION  ANALYSIS 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,L,NR,NV, IZ.MASS.NDELTA.NPSI ,NY, 
&  NRDOT.NVDOT 

DIMENSION  A(4,4) ,B(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,ALFR(4) , 
&  ALFI(4),BETA(4),WR(4),WI(4) 

C 

OPEN  (11,FILE='BIF1.RES') 

OPEN  ( 12, FILE- 'BIF2. RES') 

OPEN  (13,FILE='BIF3.RES') 


C 

C 

VehicL 

e  Parame 

IZ 

=0.0 

L 

=528 

RHO 

=  1.94 

c 

G 

=32.2 

XG 

=0.0 

MASS 

=0.0088 

U 

=3.0 

YRDOT  =  0.00000 
YVDOT  =-0.00912 
YR  =+0.00456 
YV    =-0.01434 

YPSI  =  0.01400 
YY    =  0.02000 

YDELTA=  0.00278 
NRDOT  =-0.00115 
NVDOT  =  0.00000 
NR  =-0.00296 
NV    =-0.00460 

NPSI  =  0.01000 
NY  =  -0.00250 
NDELTA=  -0.00139 

WRITE  (*,1001) 

READ  (*,*)     WNMIN,WNMAX,IWN 

WRITE  (*,1002) 

READ  (*,*)     XDMIN,XDMAX,IXD 

WRITE  0,1003) 

READ  (*,*)     ZETA 
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WRITE (*, 1100) 
READ  (*,*)     TL1 
TL1=TL1*U/L 
WRITE(*,1101) 
READ  (*,*)     TL2 
TL2=TL2*U/L 

DVR  =( IZ-NRDOT)* (MASS-YVDOT) - 
&     (MASS*XG-YRD0T) * (MASS*XG-NVD0T) 

AA11=( (IZ-NRDOT) *YPSI-(MASS*XG-YRDOT)*NPSI) /DVR 
AA12= ( ( IZ-NRDOT) *YV- (MASS*XG- YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) * (YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) *(NR-MASS*XG)) /DVR 
AA14= ( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB1  =( (IZ-NRDOT) *YDELTA- (MASS*XG- YRDOT) *NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 

DO  1  1=1, IWN 

WRITE  (*,2001)  I, IWN 
WN=WNMIN+(I-1)*(WNMAX-WNMIN)/(IWN-1) 
Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 
K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO  2  J=1,IXD 

XD=XDMIN+(J-1)*(XDMAX-XDMIN)/(IXD-1) 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 
&  AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL1,TL2) 

CALL  RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , IERR) 
DO  11  IJE=1,4 

WR(IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
11     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 
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IF  (J.GT.l)  GO  TO  10 
DEOSOO=DEOS 
XDOO  =XD 
LL=0 
GO  TO  2 
10     DEOSNN=DEOS 
XDNN  =XD 
PR=DEOSNN*DEOSOO 
IF  (PR.GT.O.DO)  GO  TO  3 
LL=LL+1 

IF  (LL.GT.3)  STOP  1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6     XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL  LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
&        AA23,AA24,BB1,BB2,A,B,TL1,TL2) 

CALL  RGG(4 ,4 , A ,B , ALFR, ALFI ,BETA ,0 , ZZZ, IERR) 
DO  12  IJE=1,4 

WR(IJE)=ALFR(IJE)/BETA(IJE) 
WI(IJE)=ALFI(IJE)/BETA(IJE) 
12     CONTINUE 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 
XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 
IF  (PRL.GT.O.DO)  GO  TO  5 
XDO=XDL 
XDN=XDM 
DEOSO=DEOSL 
DEOSN=DEOSM 
IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 
DIF=DABS (XDL-XDM) 
IF  (DIF.GT.EPS)  GO  TO  6 
XD=XDM 
GO  TO  4 
5     IF  (PRR.GT.O.DO)  STOP  3200 
XDO=XDM 
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XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF  (IL.GT.ILMAX)  STOP  3100 

DIF=DABS (XDM-XDR) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD=XDM 
4     LLL=10+LL 

WRITE  (LLL,*)  XD,WN 
3     XDOO=XDNN 

DEOSOO=DEOSNN 
2   CONTINUE 
1  CONTINUE 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless) ' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL1  (dimensional)') 

1101  FORMAT  ('  ENTER  TIME  LAG  TL2  (dimensional)') 
2001  FORMAT  (215) 

END 

SUBROUTINE  DSTABL(DEOS,WR,WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 
DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I) .LT.DEOS)  GO  TO  1 

DEOS=WR(I) 

IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS( OMEGA) 
RETURN 
END 

SUBROUTINE  LINEAR(K1 ,K2,XD,AA11 ,AA12,AA13,AA14,AA21,AA22,AA23, AA24 
&  ,BB1,BB2,A,B,TL1,TL2) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  K1,K2 

DIMENSION  A(4,4),B(4,4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 
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A(2,l 
A(2,2 

A(2,3 
A(2,4 

AC3.1 
A(3,2 
A(3,3 
A(3,4 

A(4,l 
A(4,2 
A(4,3 
A(4,4 

B(l,l 
B(l,2 
B(l,3 
B(l,4 

B(2,l 
B(2,2 
B(2,3 

B(2,4 

B(3,l 
B(3,2 
B(3,3 
B(3,4 

B(4,l 
B(4,2 

B(4,3 
B(4,4 


=AA11+BB1*K1-BB1*K1*TL2/XD-BB1*K1*TL1/XD 
=AA12-BB1*K1*TL2/XD-BB1*K1*TL1/XD 
=AA13+BB1*K2-BB1*K1*TL1+BB1*K1*TL1*TL2/XD 
=AA14+BB1*K1/XD 

=AA21+BB2*K1-BB2*K1*TL2/XD-BB2*K1*TL1/XD 
=AA22-BB2*K1*TL2/XD-BB2*K1*TL1/XD 
=AA23+BB2*K2-BB2*K1*TL1+BB2*K1*TL1*TL2/XD 
=AA24+BB2*K1/XD 

=1.0D0 
=1.0D0 
=0 . 0D0 
=0.0D0 

=1.0D0 
=0.0D0 
=0 . ODO 
=0 . ODO 

=0.0D0 

=1 . ODO- (BB1*K1*TL1*TL2/XD) 

=BB1*K2*TL1 

=0 . ODO 

=0 . ODO 

=-BB2*Kl*TLl*TL2/XD 
=1.0D0+(BB2*K2*TL1) 
=0.0D0 

=0.0D0 
=0.0D0 
=0 . ODO 
=1.0D0 


RETURN 
END 
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C     PROGRAM  HOPF.FOR 

C 

C     Hopf  Bifurcation  Analysis 

C 

C     Third  Order  Expansions:   First  Order  Approximation 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl ,K2,K3,L,NR,NV,NDELTA,NPSI ,NY,IZ,MASS, 
&  NRDOT ,  NVDOT ,  KIP ,  K2P ,  NWV ,  NVRR ,  NRVV ,  NYYY ,  NPPP 

DOUBLE  PRECISION  Mil ,M12,M13,M14,M21 ,M22,M23 ,M24, 

1  M31,M32,M33,M34,M41,M42,M43,M44, 

2  N11,N12,N13,N14,N21,N22,N23,N24, 

3  N31,N32,N33,N34,N41,N42,N43,N44, 

4  L21,L22,L23,L24,L31,L32,L33,L34, 

5  L41.L42.L43.L44 
C 

DIMENSION  A(4,4) ,T(4,4) ,TINV(4,4) ,FV1(4) , IV 1 (4) ,ZZZ(4,4) 
DIMENSION  WR(4) ,WI(4) ,TSAVE(4,4) ,TLUD(4,4) , IVLUD(4) ,SVLUD(4) 
DIMENSION  ASAVE(4,4) 


OPEN  (11,FILE='BIF1.RES') 
OPEN  ( 15, FILE=' HOPF. RES') 


c 

c 

Vehicle  Paramet< 

IZ 

=0.0 

L 

=528 

RHO 

=1.94 

G 

=32.2 

XG 

=0.0 

MASS 

=0.0088 

U 

=1.0 

c 

YRDOT 

=  0.00000 

YVDOT 

=-0.00912 

YR 

=+0.00456 

YV 

=-0.01434 

YVW 

=-0.15391 

YVRR 

=-0.05476 

YRW 

=  0.04608 

YYYY 

=  0.46800 

YPPP 

=  0.00000 

YPSI 

=  0.01400 

YY 

=  0.02000 

YDELTA=+0 . 00278 


70 


NRDOT 

=-0.00115 

NVDOT 

=   0.00000 

NR 

=-0.00296 

NV 

=-0.00460 

NVW 

=-0.00336 

NVRR 

=   0.00759 

NRVV 

=-0.05160 

NYYY 

=  0.00000 

NPPP 

=  0.00000 

NPSI 

=   0.01000 

NY 

=-0.00250 

NDELTA=-0.00139 

WRITE 

(*,1003) 

READ 

(*,*)       : 

WRITE ( 

>,1100) 

READ    ( 

:*,*) 

TL=TL*U/L 

WRITE 

(*,1006) 

READ 

(*,*)       ] 

ZETA 
TL 

DO 


DVR  = (IZ-NRDOT)* (MASS-YVDOT) - 
&     (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI ) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) * (YR-MASS) + 

&     (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

&     (MASS-YVDOT) *(NR-MASS*XG) ) /DVR 
AA14= ( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BB1  =( (IZ-NRDOT) *YDELTA-(MASS*XG-YRDOT)*NDELTA) /DVR 
BB2  = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS  =l.D-5 
ILMAX=1500 

IWN=1000 
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DO  1  11=1, IWN 

WRITE  (*,2001)  II 

READ(11,*)  XD,WN 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 

K3=K2 

C       Start  Hopf  Bifurcation  Analysis 

C 

C       Evaluate  Nonlinear  Rudder  Expansion  Coefficients 

C 

K2=0 . 0D0 

K1P=K1-K1*TL*U/XD 
K2P=K2-K1*TL/XD 
C 

DPPV=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K2P 

&       +  0.5D0*K1*TL/XD  +  3.D0*K1*(TL*U) **3/(3 .D0*XD**3) 
DPVV=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K2P*K2P 

&       +  3.D0*U*TL*TL*TL*K1/(3.D0*XD**3) 

DPPR=-  (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K3 
DPRR=-  (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K3*K3 
DPPY=-  (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K1/XD 

&       -  3.D0*TL*TL*U*U*K1/(3.D0*XD**3) 

DPYY=-  (1 .DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1*K1/ (XD**2) 

&       +  3.D0*TL*U*K1/(3.D0*XD**3) 

DVVR=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K2P*K2P*K3 
DVRR=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K2P*K3*K3 
DVVY=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . DO*K 1*K2P*K2P/XD 

&       -  3.D0*K1*TL*TL/(3.D0*XD**3) 

DVYY=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K1*K2P/ (XD**2) 

&       +  3.D0*TL*K1/(3.D0*XD**3) 

DRRY=-  ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K3*K3/XD 
DRYY=-  (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K1*K3/ (XD**2) 
DPVR=-  (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K2P*K3 
DPVY=-  (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K1*K2P/XD 

&       -  6.D0*TL*TL*U*K1/(3.D0*XD**3) 

DPRY=-  ( 1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K1*K3/XD 
DVRY=-  ( 1 . DO/ (3 . D0*D0**2) ) *6 . DO*Kl *K2P*K3/XD 
DPPP=-  ( 1 . DO/ (3 . D0*D0**2) ) *1 . D0*K1P*K1P*K1P 

&       +  K1*TL*U/(6.D0*XD)  +  (K1*(TL*U)**3)/(3.D0*XD**3) 
DVVV=-  ( 1 . DO/ (3 . D0*D0**2) ) *1 . D0*K2P*K2P*K2P 

&       +  K1*(TL**3)/(3.D0*XD**3) 

DRRR=-  ( 1 . DO/ (3 . D0*D0**2) ) *1 . D0*K3*K3*K3 

DYYY=-  (1.D0/(3.D0*D0**2))*(K1/XD)**3-K1/(3.D0*XD**3) 

&       -  K1/(3.D0*XD**3) 
C 

C       Evaluate  Transformation  Matrix  of  Eigenvectors 
C 
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CALL  LINEAR (Kl ,K3,XD,AA11 ,AA12,AA13,AA14,AA21 ,AA22,AA23 , 
&    AA24,BB1,BB2,A,TL) 

DO  11  1=1,4 
DO  12  J=l,4 

ASAVE(I,J)=A(I,J) 

12  CONTINUE 
11   CONTINUE 

CALL  RG(4,4,A,WR,WI,1,ZZZ,IV1,FV1,IERR) 
CALL  DS  OMEG ( IEV , WR , WI , OMEGA , CHECK ) 
OMEGAO=OMEGA 
DO  50  1=1,4 

T(I,1)=  ZZZ(I.IEV) 
T(I,2)=-ZZZ(I,IEV+1) 
50   CONTINUE 

IF  (IEV.EQ.l)  GO  TO  13 
IF  (IEV.EQ.2)  GO  TO  14 
IF  (IEV.EQ.3)  GO  TO  15 
STOP  3004 

14  DO  60  1=1,4 

T(I,3)=ZZZ(I,1) 
T(I,4)=ZZZ(I,4) 
60   CONTINUE 
GO  TO  17 

15  DO  70  1=1,4 

T(I,3)=ZZZ(I,1) 
T(I,4)=ZZZ(I,2) 
70   CONTINUE 
GO  TO  17 

13  DO  16  1=1,4 

T(I,3)=ZZZ(I,3) 
T(I,4)=ZZZ(I,4) 

16  CONTINUE 

17  CONTINUE 
C 

C       Normalization  of  the  Critical  Eigenvector 
C 

CALL  NORMAL (T) 
C 

C       Definition  of  Mij 
C 

M11=T(1,1) 

M21=T(2,1) 

M31=T(3,1) 

M41=T(4,1) 

M12=T(1,2) 

M22=T(2,2) 
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M32=T(3,2) 

M42=T(4,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 

M43=T(4,3) 

M14=T(1,4) 

M24=T(2,4) 

M34=T(3,4) 

M44=T(4,4) 

c 

c 

Definition  of  Lij 

c 

L21=  DPPV*M11*M11*M21 

& 

+  DPRR*M11*M31*M31 

& 

+  DVVR*M21*M21*M31 

& 

+  DVYY*M21*M41*M41 

& 

+  DPVR*M11*M21*M31 

& 

+  DVRY*M21*M41*M31 

& 

+  DRRR*M31*M31*M31 

+  DPVV*M11*M21*M21  +  DPPR*M11*M11*M31 

+  DPPY*M11*M11*M41  +  DPYY*M11*M41*M41 

+  DVRR*M21*M31*M31  +  DVVY*M21*M21*M41 

+  DRRY*M31*M31*M41  +  DRYY*M31*M41*M41 

+  DPVY*M11*M21*M41  +  DPRY*M11*M41*M31 

+  DPPP*M11*M11*M11  +  DVVV*M21*M21*M21 
+  DYYY*M41*M41*M41 


PPV  =  M11*M11*M22  +  2.0*M11*M12*M21 
PVV  =  M12*M21*M21   +  2.0*M11*M21*M22 
PPR  =  M11*M11*M32  +  2.0*M11*M12*M31 
PRR  =  M12*M31*M31   +   2.0*M11*M31*M32 
PPY  =  M11*M11*M42  +  2.0*M11*M12*M41 
PYY  =  M41*M41*M12  +  2.0*M11*M41*M42 
WR  =  M21*M21*M32  +   2.0*M31*M21*M22 
VRR  =  M22*M31*M31   +  2.0*M31*M32*M21 
VVY  =  M21*M21*M42  +  2. 0*M41*M21*M22 
VYY  =  M22*M41*M41   +  2. 0*M41*M42*M21 
RRY  =  M31*M31*M42  +  2. 0*M41*M31*M32 
RYY  =  M32*M41*M41   +  2. 0*M41*M42*M31 
PVR  =  M11*M21*M32+M31*(M11*M22+M12*M21) 
PVY  =  M11*M21*M42+M41*(M11*M22+M12*M21) 
PRY  =  M11*M41*M32+M31*(M11*M42+M12*M41) 
VRY  =  M21*M41*M32+M31*(M21*M42+M22*M41) 
PPP  =  3.0*M11*M11*M12 
VVV  =  3.0*M21*M21*M22 
RRR  =  3.0*M31*M31*M32 
YYY  =  3.0*M41*M41*M42 

L22=DPPV*PPV+DPVV*PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*VVV 
+DRRR*RRR+DYYY*YYY 

PPV  =  M12*M12*M21   +  2.0*M11*M12*M22 
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PVV  =  M11*M22*M22  +  2. 0*M12*M21*M22 
PPR  =  M12*M12*M31   +  2.0*M11*M12*M32 
PRR  =  M11*M32*M32  +  2.0*M12*M31*M32 
PPY  =  M12*M12*M41   +  2.0*M11*M12*M42 
PYY  =  M11*M42*M42  +  2. 0*M12*M41*M42 
WR  =  M22*M22*M31   +  2. 0*M21*M22*M32 
VRR  =  M21*M32*M32  +   2. 0*M22*M31*M32 
VVY  =  M22*M22*M41   +  2. 0*M21*M22*M42 
VYY  =  M21*M42*M42   +  2. 0*M22*M41*M42 
RRY  =  M32*M32*M41   +  2. 0*M31*M32*M42 
RYY  =  M31*M42*M42  +  2.0*M32*M41*M42 
PVR  =  M12*M22*M31   +  M32*(M11*M22+M12*M21) 
PVY  =  M12*M22*M41   +  M42* (M11*M22+M12*M21) 
PRY  =  M12*M42*M31   +  M32*(M11*M42+M12*M41) 
VRY  =  M22*M42*M31   +  M32*(M21*M42+M22*M41) 
PPP  =  3.0*M11*M12*M12 
WV  =  3.0*M21*M22*M22 
RRR  =  3.0*M31*M32*M32 
YYY  =  3.0*M41*M42*M42 

L23=DPPV*PPV+DPVV*PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
&  +DWR*VVR+DVRR*VRR+DVVY*WY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

&  +DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*VVV 

&  +DRRR*RRR+DYYY*YYY 

L24=  DPPV*M12*M12*M22  +  DPVV*M12*M22*M22  +  DPPR*M12*M12*M32 

&  +  DPRR*M12*M32*M32  +  DPPY*M12*M12*M42  +  DPYY*M12*M42*M42 

&  +  DVVR*M22*M22*M32  +  DVRR*M22*M32*M32  +  DVVY*M22*M22*M42 

&  +  DVYY*M22*M42*M42  +  DRRY*M32*M32*M42  +  DRYY*M32*M42*M42 

&  +  DPVR*M12*M22*M32  +  DPVY*M12*M22*M42  +  DPRY*M12*M42*M32 

&  +  DVRY*M22*M42*M32  +  DPPP*M12*M12*M12  +  DVVV*M22*M22*M22 

&  +  DRRR*M32*M32*M32   +  DYYY*M42*M42*M42 

L31=L21 
L32=L22 
L33=L23 
L34=L24 

L21=L21*BB1*U*U 

L22=L22*BB1*U*U 

L23=L23*BB1*U*U 

L24=L24*BB1*U*U 

L31=L31*BB2*U*U 

L32=L32*BB2*U*U 

L33=L33*BB2*U*U 

L34=L34*BB2*U*U 

L41=-0 . 5*M11*M11* (M21+U*Mll/3 . 0) 

L42=-M11* (M12*M21+0 . 5*Mll*M22+0 . 5*U*M12*M11) 
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L43=-M12* (Ml l*M22+0 . 5*M12*M21+0 . 5*U*M1 1*M12) 
L44=-0 . 5*M12*M12* (M22+U*M12/3 . 0) 

C11=(1/6.0)*YVW*(M21*M21*M21)+0.5*YVRR*(M31*M31*M21)  + 
&  0 . 5*YRVV*M21*M21*M31+ ( 1/6 . 0) *YPPP*Mll*Mll*Mll+(l/6 . 0) *YYYY*M41*M41*M41 

C12=  (1/6 .  0)  *YVW*  (3*M21*M21*M22)+0 .  5*YVRR*  (M31*M31* 
&  M22+2*M31*M32*M21)+0.5*YRVV*(M21*M21*M32+2*M31*M21*M22)+(1/6.0)*YPPP*3*M11*M11*M12 
&  +(1/6.0)*YYYY*3*M41*M41*M42 


C13=  ( 1/6 .  0)  *YVW*  (3*M21*M22*M22)+0 .  5*YVRR*  (M21*M32* 
&  M32+2*M31*M32*M22)+0.5*YRVV*(M22*M22*M31+2*M32*M21*M22)+(1/6.0)*YPPP*3*M11*M12*M12 
&  +(1/6.0)*YYYY*3*M41*M42*M42 

C14=  (1/6.0)  *YVW*  (M22*M22*M22)  +0 .  5*YVRR*  (M32*M32*M22)  + 
&  0 . 5*YRVV*M22*M22*M32+ (1/6.0) *YPPP*M12*M12*M12+ (1/6.0) *YYYY*M42*M42*M42 

C21=(1/6.0)*NVVV*(M21*M21*M21)+0.5*(NVRR*M31*M31*M21)+ 
&   0.5*NRVV*M21*M21*M31+(1/6.0)*NPPP*M11*M11*M11+(1/6.0)*NYYY*M41*M41*M41 

C22=  ( 1/6 . 0)  *NVW*  (3*M21*M21*M22)+0 .  5*NVRR*  (M31*M31* 
&  M22+2*M31*M32*M21)+0 . 5*NRVV* (M21*M21*M32+2*M31*M21*M22) + (1/6 . 0) *NPPP*3*M11*M11*M12 
&  +(1/6.0)*NYYY*3*M41*M41*M42 

C23=  (1/6 . 0)  *NVW*  (3*M21*M22*M22)+0  .  5*NVRR*  (M21*M32* 
&  M32+2*M31*M32*M22)+0.5*NRVV*(M22*M22*M31+2*M32*M21*M22)+(1/6.0)*NPPP*3*M11*M12*M12 
&  +(1/6.0)*NYYY*3*M41*M42*M42 

C24=  (1/6.0)  *NVW*  (M22*M22*M22)  +0 .  5*NVRR*M32*M32*M22+ 
&   0 . 5*NRVV*M22*M22*M32+ (1/6.0) *NPPP*M12*M12*M12+ (1/6.0) *NYYY*M42*M42*M42 


D11=(C11*(IZ-NRD0T)+C21*(YRD0T-M*XG))/DVR 
D12= (C12* (IZ-NRDOT) +C22* (YRDOT-M*XG) ) /DVR 
D13= (C13* (IZ-NRDOT) +C23* (YRDOT-M*XG) ) /DVR 
D14= (C14* (IZ-NRDOT) +C24* (YRDOT-M*XG) ) /DVR 

D21= (CI 1* (NVDOT-M*XG) +C21* (M-YVDOT) ) /DVR 
D22= (C12* (NVDOT-M*XG) +C22* (M-YVDOT) ) /DVR 
D23= (C13* (NVDOT-M*XG) +C23* (M-YVDOT) ) /DVR 
D24= (C14* (NVDOT-M*XG) +C24* (M-YVDOT) ) /DVR 

C 

C       Invert  Transformation  Matrix 

C 

DO  20  1=1,4 
DO  30  J=l,4 


76 


TINV(I,J)=0.0 
TSAVE(I,J)=T(I,J) 
30     CONTINUE 
20   CONTINUE 

CALL  DLUD(4,4,TSAVE,4,TLUD,IVLUD) 
DO  40  1=1,4 

IF  (IVLUD(I).EQ.O)  STOP  3003 
40   CONTINUE 

CALL  DILU(4,4,TLUD,IVLUD,SVLUD) 
DO  80  1=1,4 
DO  90  J=l,4 

TINV(I,J)=TLUD(I,J) 
90     CONTINUE 
80   CONTINUE 
C 

C       Check  Inv(T)*A*T 
C 

IMULT=2 

IF  (IMULT.EQ.l)  CALL  MULT(TINV,ASAVE,T) 
IF  (IMULT.EQ.O)  STOP 
C 

C       Definition  of  Nij 
C 

N11=TINV(1,1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
N41=TINV(4,1) 
N12=TINV(1,2) 
N22=TINV(2,2) 
N32=TINV(3,2) 
N42=TINV(4,2) 
N13=TINV(1,3) 
N23=TINV(2,3) 
N33=TINV(3,3) 
N43=TINV(4,3) 
N14=TINV(1,4) 
N24=TINV(2,4) 
N34=TINV(3,4) 
N44=TINV(4,4) 
C 

R11=N12*L21+N13*L31+N14*L41+N12*D11+N13*D21 
R12=N12*L22+N13*L32+N14*L42+N12*D12+N13*D22 
R13=N12*L23+N13*L33+N14*L43+N12*D13+N13*D23 
R14=N12*L24+N13*L34+N14*L44+N12*D14+N13*D24 
R21=N22*L21+N23*L31+N24*L41+N22*D11+N23*D21 
R22=N22*L22+N23*L32+N24*L42+N22*D12+N23*D22 
R23=N22*L23+N23*L33+N24*L43+N22*D13+N23*D23 
R24=N22*L24+N23*L34+N24*L44+N22*D14+N23*D24 
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c 

C       Evaluate  Alpha'  and  Omega' 
C 

DINC=0.001 

XDR  =XD+DINC 

XDL  =XD-DINC 

XD  =XDR 

CALL  LINEAR(K1,K3,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
&    AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
ALPHR=DEOS 
OMEGR=FREQ 

XD=XDL 
C 

CALL  LINEAR (K1,K3,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
&    AA24,BB1,BB2,A,TL) 

CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 
C 

DALPHA= (ALPHR-ALPHL) / (XDR-XDL) 

DOMEGA= (OMEGR-OMEGL) / (XDR-XDL) 
C 

C       Evaluation  of  Hopf  Bifurcation  Coefficients 
C 

C0EF1=3 . 0*R1 1+R13+R22+3 . 0*R24 

C0EF2=3 . 0*R2 1+R23-R12-3 . 0*R14 

AMU2  =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*C0EF1 

TAU2  =- (C0EF2-D0MEGA*C0EF1/DALPHA) / (8 . 0*0MEGA0) 

PER     =2.0*3. 1415927/0MEGA0 

PER     =PER*U/L 
C 

WRITE  (15,2002)  XD , WN , COEF 1 , DALPHA , PER 


1  CONTINUE 


1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  WN  (dimensionless) ' ) 

1002  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  XD  (dimensionless)') 

1003  FORMAT  ('  ENTER  DAMPING  RATIO') 
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1006  FORMAT  ('  ENTER  DSAT  ') 

1100  FORMAT  ('  ENTER  TIME  LAG  TL  (dimensional) ') 

2001  FORMAT  (215) 

2002  FORMAT  (5E15.5) 


END 


SUBROUTINE  DSTABL(DEOS,WR,WI , OMEGA) 
C 

C     EVALUATES  THE  EIGENVALUE  WITH  THE  MAXIMUM  REAL  PART 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR(4),WI(4) 
DE0S=-1.0D+20 
DO  1  1=1,4 

IF  (WR(I).LT.DEOS)  GO  TO  1 
DEOS=WR(I) 
IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 
C 

SUBROUTINE  LINEAR(K1,K2,XD,AA11 ,AA12,AA13,AA14,AA21, 
&  AA22,AA23,AA24,BB1,BB2,A,TL) 
C 
C     FORMS  THE  LINEARIZED  MATRIX  A  (time  delay  1st  order  approximation 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  K1,K2 

DIMENSION  A(4,4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2, 1)=AA11+BB1*K1-BB1*K1*TL/XD 

A(2,2)=AA12-BB1*K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A(2,4)=AA14+BB1*K1/XD 

A(3 , 1) =AA21+BB2*K1-BB2*K1*TL/XD 

A(3 , 2) =AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A(3,4)=AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 
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A(4,3)=0.0D0 
A(4,4)=0.0D0 
RETURN 
END 


SUBROUTINE  DSOMEG ( I JK , WR , WI , OMEGA , CHECK) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR(4),WI(4) 
CHECK"- i.OD+25 
DO  1  1=1,4 

IF  (WR(I) .LT. CHECK)  GO  TO  1 

CHECK=WR(I) 

IJ=I 
CONTINUE 

OMEGA=DABS(WI(IJ)) 
IF  (WI(IJ) .GT.O.DO)  IJK=IJ 
IF  (WI(IJ) .LT.O.DO)  IJK=IJ-1 
RETURN 
END 


SUBROUTINE  NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  T(4,4) ,TNOR(4,4) 

CFAC=T(1,1)**2+T(1,2)**2 

IF  (DABS(CFAC).LE.(l.D-lO))  STOP  4001 

TN0R(1,1)=1.D0 

TN0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 

TN0R(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 

TN0R(4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2))/CFAC 

TN0R(1,2)=0.D0 

TN0R(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 

TN0R(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 

TN0R(4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2))/CFAC 

DO  1  1=1,4 
DO  2  J=l,2 

T(I,J)=TNOR(I,J) 
2   CONTINUE 
1  CONTINUE 

RETURN 
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END 
c======================================================: 

SUBROUTINE  MULT(TINV,A,T) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  TINV(4,4) ,A(4,4) ,T(4,4) ,A1(4,4) ,A2(4,4) 
DO  1  1=1,4 
DO  2  J=l,4 

Al(I,J)=O.DO 

A2(I,J)=0.D0 

2  CONTINUE 
1  CONTINUE 

DO  3  1=1,4 
DO  4  J=l,4 
DO  5  K=l,4 

A1(I,J)=A(I,K)*T(K,J)+A1(I,J) 

5  CONTINUE 
4   CONTINUE 

3  CONTINUE 
DO  6  1=1,4 

DO  7  J=l,4 
DO  8  K=l,4 

A2(I,J)=TINV(I,K)*A1(K,J)+A2(I,J) 
8     CONTINUE 
7   CONTINUE 

6  CONTINUE 

DO  11  1=1,4 
C        WRITE  (*,101)  (A(I,J) ,J=1,4) 

11  CONTINUE 

DO  12  1=1,4 
C        WRITE  (*,101)  (T(I,J),J=1,4) 

12  CONTINUE 

DO  10  1=1,4 

WRITE  (*,101)  (A2(I,J),J=1,4) 
10  CONTINUE 
C      WRITE  (*,101)  A2(l,l) 
RETURN 
101  FORMAT  (4E15.5) 
END 


81 


82 


LIST  OF  REFERENCES 


1.  Beck,  R.  F.  [1976]  "Forces  and  moments  on  a  ship  moving  in  a  canal", 
Report  No.  179,  Dept.  of  Naval  Architecture  and  Marine  Engineering, 
The  University  of  Michigan,  Ann  Arbor. 

2.  Bernitsas,  M.  M.  and  Kekridis,  N.  S.  [1984]  "Nonlinear  simulation  of  time 
dependent  towing  of  ocean  vehicles",  Report  No.  283,  Dept.  of  Naval 
Architecture  and  Marine  Engineering,  The  University  of  Michigan,  Ann 
Arbor. 

3.  Chow,  S.-N.  and  Mallet-Paret,  J.  [1977]  "Integral  averaging  and  bifur- 
cation", Journal  of  Differential  Equations,  26,  pp.  112-159. 

4.  Clayton,  B.  R.  and  Bishop,-  R.  E.  D.  [1982]  Mechanics  of  Marine  Vehicles 
(Gulf  Publishing  Company,  Houston). 

5.  Comstock,  J.  P.  [1977]  Principles  of  Naval  Architecture  (The  Society  of 
Naval  Architects  and  Marine  Engineers,  New  York). 

6.  Guckenheimer,  J.  and  Holmes,  P.  [1983]  Nonlinear  Oscillations,  Dynam- 
ical Systems,  and  Bifurcations  of  Vector  Fields  (Springer- Verlag,  New 
York). 

7.  Hassard,  B.  and  Wan,  Y.H.  [1978]  "Bifurcation  formulae  derived  from 
center  manifold  theory" ,  Journal  of  Mathematical  Analysis  and  Applica- 
tions, 63,  pp.  297-312. 

8.  Venne,  D.  V.  [1992]  "Effects  of  positional  information  time  lags  on  motion 
stability  of  autonomous  vehicles" ,  Master's  Thesis,  Naval  Postgraduate 
School,  Monterey. 


83 


84 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 

1.  Defense  Technical  Information  Center  2 

8725  John  J.  Kingman  Rd.  STE  0944 

Ft.  Belvoir,  VA  22060-6218 

2.  Dudley  Knox  Library  2 

Naval  Postgraduate  School 

411  Dyer  Rd. 

Monterey,  California  93943-5101 

3.  Chairman,  Code  ME 1 

Department  of  Mechanical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93943 

4.  Professor  Fotis  A.  Papoulias,  Code  ME/Pa  3 

Department  of  Mechanical  Engineering 

Naval  Postgraduate  School 
Monterey,  California  93943 

5.  Panos  E.  Kapasakis 3 

Argyriou  Mathesi  4 

Salamina  18900 
Greece 

6.  Embassy  of  Greece 1 

Naval  Attache 

2228  Massachusetts  Avenue,  N.W. 
Washington,  D.C.  20008 

7.  Naval  Engineering  Curricular  Office,  Code  34 1 

Naval  Postgraduate  School 

Monterey,  California  93943 


85 


15     483NPG 
TH 

10/99  22527-200 


MOW.  'eRAOU  ^..^ 


.1 


